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ABSTRACT 

The basic physical relationships involved in control of a flexible air- 
craft disturbed by random wind gusts were used in formulating the sur- 
face location problem as one in optimal control of a distributed system, 
using a limited number of point-force controllers. The three phases 
of this problem --estimation, control, and surface placement- -were 
then solved by means of the matrix minimum principle and the calculus 
of variations. The variational equations were greatly simplified and 
the order of the problem considerably reduced by the use of optimal 
controllers at each stage in the search for optimal surface locations. 
These simplifications, plus advanced computational techniques made 
the general solution of the problem practically feasible. 

Aircraft physics had to be investigated in great detail in order to ob- 
tain general equations expres sing the distributed nature of the system 
exactly. A computer program was written which stored these equations 
and used them in sovling the surface location problem for a general 
aircraft’. 

This program was tested on a fourteenth order model of the Lockheed 
C-5A transport aircraft. As a guide for future applications, the deri- 
vation of the model parameters was carried through explicitly. The 
results of the optimization study were then analyzed in an attempt to 
recognize and develop a "strategy 1 ' for control surface placement. 

A simple practical strategy was developed for systems with stress and 
stress rate responses and slightly simplified physics. This strategy, 
which was also verified for the C-5A model, has the two major ad- 
vantages of (1) insight into the tradeoffs involved in surface location, 
and (2) partial insight into global solutions. 

The major contributions of this thesis are (1) a computationally feasible 
general solution to the surface placement problem, (2) a practical ap- 
plication of optimal flexure control to a trial model of the C-5A trans- 
port, (3) implementation of efficient computational techniques for 
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solving state covariance and Riccati equations, and (4) recognition of 
the nature of the optimal solutions and presentation of a search 
strategy which makes use of the basic aircraft flexure physics. 
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INTRODUCTION 


Historically the task of designing a flight control system has not 
included choosing the sizes or locations of the aerodynamic control 
surfaces. In many flight control problems the control designer has not 
even been free to choose the actuators which drive these surfaces. 
There has not been a need to give the designer these freedoms; sur- 
face locations and control authority are chosen for rigid body trim 
and maneuverability, and flight control designers have been designing 
rigid body control systems . 

The advent of active flexure control has changed this picture . It 
is obvious from the flexure mode physics that no mode can be effec- 
tively controlled through a surface located near a node point. Outboard 
ailerons, for example, are evidently more effective for control of a 
wing torsion mode than elevator flaps. Clearly each directly con- 
trolled mode must be controllable, and some force producer locations 
will be more effective than others . 

The freedom to choose control surface locations creates several 
Tiew issues: How many control surfaces should be used; what are the 
trade-offs between them; can physically intuitive guidelines for sur- 
face placement be found; how much performance will added surfaces 
purchase, or will optimal positions change radically when new sur- 
faces are added? These issues can be reduced to two fundamental 
questions : 

1 . What intuitive considerations (if any) are most significant 

in choosing control surface locations -- what are the funda- 
mentals ? 
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2. What degree of over-all performance improvement can be 
achieved by using the freedom to choose surface locations? 

The primary concern of this thesis will be to answer the first question- - 
how to identify the fundamentals of control surface location. By 
closely examining the flexure physics of the aircraft, one can under- 
stand the constraints posed by a limited number of control surfaces 
for a large number of vibration modes and develop a strategy for 
roughly estimating optimal locations. This rough estimation may be 
seen either as method of choosing initial locations for a search pro- 
cedure, or as a "rule of thumb" to be used in lieu of an exact solution 
to the problem. In both cases the precise optimization problem serves 
as a check on the validity of the approximations used. 

The C-5A transport (see Fig. 1) is an excellent vehicle for such 
investigations; structural flexure modes, being a significant factor in 
performance evaluation, can cause some control surface locations to 
be definitely superior to others. The great size of the C-5A makes it a 
ready candidate for multiple surfaces of many types -- ailerons, 
spoilers, canards, leading edge flaps, elevators, rudders, etc. Be- 
cause of the many modes (15 modes have been calculated for both 
pitch and lateral models) and the possibility of so many types of multi- 
ple surfaces, the C-5A surface location problem definitely demands a 
control -theoretic treatment . The scale of the problem also permits 
proposed fundamentals to be adequately tested. 

Optimal control theory will be shown to provide a qualitative 
framework for analyzing the performance of an aircraft with any speci- 
fied set of control surfaces and surface locations. The success of 
programs such as Honeywell's Load Alleviation and Mode Stabilization 
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Fig. 1 Lockheed C-5A Transport Aircraft 
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System for the B-52 attests to the practicality of quadratic criteria 
in evaluating performance. Hence, the second question above may be 
answered by simply casting the relevant equations of the aircraft into 
an optimal control framework and evaluating the system performance 
by such measures as stresses and mean square accelerations. 

The body of this thesis is organized around the three contri- 
butions which it makes: 

1 . Solution of the aerodynamic surface placement problem 

as a distributed- system problem in optimal control theory. 

2. Practical implementation of the solution procedure; this in- 
cludes a new representation of control forces as a function of 
position and description of a fast state-of-the-art computer 
program applicable to a wide variety of aircraft. 

3. A strategy of surface placement -- an approximate method of 
choos ing locations for a restricted class of problems. 

The first section motivates the modelling process by describing 
the fundamentals of the physical system involved. Section II gives the 
formulation and solution of the theoretical problem. Section III de- 
scribes implementation of the technique -- a precise mathematical 
model of the physics of typical aircraft, and a computational scheme 
for carrying out the desired calculations. Application of the computer 
program to a scaled-down model of a large transport aircraft, the 
C-5A, is given in Section IV. Section V explains a possible strategy 
of surface placement applicable to an important but restricted class of 
problems. The sections are written with a maximal amount of in- 
dependence in order to accommodate readers with varying interests . 



SECTION I 


THE PHYSICAL SYSTEM 

The introduction discusses surface location in an intuitive fashion-- 
indeed the existence and nature of the control surface location problem 
are evident on such grounds. It is not evident, however, how to ap- 
proach and organize the problem- -in what ways does the solution de- 
pend on pure flexure properties of the aircraft, on dynamic properties of 
the rigid vehicle, on the probabilistic structure of gusts which strike the 
plane, or on the responses one chooses to evaluate . This section 
makes clear the physical relationships between these variables and 
motivates a state -variable description of the problem. 

A. BASIC SYSTEM STRUCTURE 

Because the aerodynamic forces which govern the motion of an 
aircraft in flight are basically nonlinear (Bernoulli^ principle pre- 
dicts a square-law dependence of forces on relative wind velocity), a 
reasonable control problem must treat the vehicle in one or more . 
specified M flight conditions.' 1 A flight condition is generally specified 
by a steady velocity of the aircraft and its orientation with respect to 
the mean wind striking the vehicle, as well as the settings of throttle 
and control surfaces which maintain this condition of steady flight. 

The remainder of the thesis presumes that the primary n plant ,T is a 
large flexible aircraft whose motions are perturbed about such a 
steady flight condition. This assumption is reasonable since such 
vehicles are designed for long periods of steady flight. 

The perturbing forces acting on an aircraft in steady flight are 
(1) random gust disturbances (relative to the mean wind, about which/ 
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all aerodynamic forces are linearized), and (2) control forces resulting 
from deflections of the aerodynamic surfaces* Because of its in- 
efficiency, engine control is not considered in detail, although it can 
be readily incorporated into the framework of analysis presented. 

The random excitations of the plant due to the gusts give rise to 
several undesirable responses, such as stresses and distributed torques 
which cause structural fatigue of the wings and tail, and accelerations 
affecting ride quality and safety of cargo. The objective of flexure 
control is to produce control forces in such a way as to counteract 
these undesirable effects of the stochastic gust inputs* Generally the 
control engineer is aware of certain locations on the vehicle where 
large responses are least desirable, e.g. , weak points in the basic 
structure of the aircraft, crucial points such as those where engines 
are fixed to the wings, or passenger stations in the fuselage. In 
light of this information, response locations are assumed to be pre- 
specified independently of control surface locations. This is possible 
because the basic structural properties of the vehicle are independent 

of surface location; aircraft designers specifically avoid major 

» 

9 

structural' modifications for control surfaces, as these tend to weaken 
the overall vehicle* (Contrary to intuitive reasoning, it will be shown 
in Section III that the incidence of control forces is not at the hinge 
line of a flap --hence there are even more basic physical reasons for 
choosing response locations independently of control surface locations). 

A fundamental though unproven tenet of control theory is that in 
order to control something, it must be measured. The problem of 
precisely what signals to sense appears to be quite complicated, as it 
involves a coupling of control and estimation problems which yields 
computationally intractable mathematics (with present technology), 
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so this thesis henceforth assumes sensor locations to be prespecified. 
Possible types of sensors include rate and position gyros, accelero- 
meters, angle of attack probes, etc. These sensors pick up certain 
combinations of dynamic state and control variables such as flexure 
mode accelerations, pitch rate, control forces generated by aileron 
displacement, etc. An important distinction should be made between 
responses and sensed signals: sensor locations maybe entirely dif- 

ferent than response locations --in fact there is quite possibly no 
* * * * 

necessary connection whatsoever between responses and measured 
signals. Control of responses depends on their functional relation to 
the state of the system, and not on measurements of their exact values 
at a given time. Due to modelling uncertainty and component vari- 
ations, sensed signals are customarily assumed to be corrupted to 
some extent. Sensor dynamics are generally ignored, as they are 
several orders of magnitude faster than aerodynamic and structural 
responses . 

Measured signals, which often give only a partial measure of 
the state of the aircraft, are all that the control designer may use in 
constructing feedback signals to the control surfaces. Furthermore, 
hardware constraints demand that the control law have a minimum of 
dynamics and a relatively simple form. In effect this dictates (at 
most) constant gains on feedback of system states. The existence of 
adequate (though sub optimal) controllers with constant element values 
and relatively simple form convinces one that optimal controllers 
should be able to perform effectively within the same constraints. 

The freedom to choose surface locations essentially expands the 
range of forces and moments that the control designer may exert on 
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the vehicle in order to damp oat its responses. Presumably, this new 
freedom should allow him to do a better job with a specified set of 
responses and measured signals . 

Figure 1 schematizes the basic system structure as described 
above . 

B. MATHEMATICAL MODELLING OF THE SYSTEM 

The most fundamental assumption underlying the modelling pro- 
cess is linearity. Linearity, in view of the valuable body of infor- 
mation on linear optimal control, is certainly desirable on a mathe- 
matical basis ; the definition of the problem in the framework of a 
given flight condition has also made linearity a justifiable physical 
assumption. It should be noted that linearization has its costs-- 
practical system designs require examination of many flight con- 
ditions and great expense in deriving coefficients for many linear 
models. If one makes the claim, however, that fatigue lifetimes are 
governed primarily by the constant flexing of an aircraft in normal 
(gusty) flight, linearized analyses are also a good way to attack the 
flexure control problem. In view of the fact that aerodynamic stresses 
over the cross section of a wing may vary by several tons in an average 
wind gust, this claim seems quite reasonable. 

Another major assumption is that the distributed nature of the 
system may be adequately approximated by the use of a finite (and 
reasonably small) number of flexure modes. Although this assumption 
is well -justified in this problem (due to aerodynamic damping of higher 
modes), it poses formidable theoretical problems for distributed 
systems in general. 

On the basis of the foregoing physical description and the as- 
sumption of linearity,, one may model the functional relationships of a 
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complete system by the equations and variables below. 

Definition of Variables 

x^ rigid body state {angle of attack, sideslip, roll, 
pitch and yaw rates, etc.) 
x ^ elastic state (mode deflections and rates) 

5-3 gust velocities generated by gust filter 
x 4 gust forces (aerodynamic effect of gusts) 

££5 control forces (aerodynamic effect of control surface 
deflections --may be taken to include actuator dy- 
namics) 

c control signals (actuator drives) 

m sensor measurements 

_r responses to be minimized 

white noise inputs to gust filter 
B 2 wh i te noise inputs to sensor measurement s 
Y. vector of control surface locations 

System Equations 

+ — 1—2 + —1—4 + 5 (Rigid Body Equations 

of Motion) (1) 

x 2 = L 2 2£j + ^2—2 + —2—4 + —2^—5 (El exure Equations) (2) 

S 3 = (Gust Filter) (3) 

X4 = F 4 X 4 + S 3 (Gust Aerodynamics) (4) 

x 5 = F^x,- + c (Control Surface Aerodynamics) (5) 

m - AjX j + 

v_ = I-I j[ x 1 -f H 2 x 2 4 H4X4 + Dr-(y)x^ (Response Equations) (7) 
= + —z—2 + —4—4 + — 5^—5 -5 + P 5 (X)£ 


4a 


+ A .x. + 
-4-4 


-2-^2 (Censor Equations) 


(6) 
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Sta te a ugmentat ion 


Define the new augmented states: 



The vector x is taken to be the optimal linear estimator 
— o 

of the vector x q based on the measured signals m (to be 
explained later) , With these definitions, Eqs * 1 through 7 


become : 





m = A x + B ri 
r = Hx + D c 




A = [ A x A 2 0 A 4 0] B = [0 B 2 ] 

H = [Hj H 2 0 H 4 D 5 (j)F 5 ] D = [D 5 (y)] 


Control Problem 

According to the model, one has perfect knowledge of the state 
x,. {aerodynamic forces due to control signals), since , the control 
signals, are generated directly by whatever feedback system is used. 
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These signals are the only inputs to So given c, one can in 

theory calculate x and if any measurements were to contain x^, 
it could be subtracted out exactly to obtain a set of modified signals 
independent of x ^ ; note that this assumption has been implemented in 
Eq, 6 by making m independent of x^. 

If one knew the remaining states of the system, w °uld 

wish to use a control law of the form 

£ = K(z)i (U) 

since x by definition contains all information on the state of the sys- 
tem, K(y;) is chosen to be constant rather than time -varying because 
a time -varying gain would be too difficult to implement; it will later 
be shown that for long flight times (which is certainly the case), a 
constant gain K(y) is optimal at any rate. 

But x^ is not completely known, in general, since the system 
can only be measured by a limited number of stochastically disturbed 
sensors, for economic and practical reasons. Therefore, one hopes 
to use a linear estimat e , x , which (in a way to be defined in Section II) 
best approximates x^. Intuition readily tells one that must 

satisfy an equation of the form 

X =Fx+Gm* 

— o * O “O — o 

since this is the most general linear dynamic system making use of 
the measurements available. In practice, this system (the Kalman 
filter) will not be easy to build. It will be used in this thesis, however, 
because of its optimality properties ; it provides a justifiable standard 

If m should include x the author has shown that the Kalman 
filter will perfectly estimate this state anyway. 
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by which to compare end results of the location problem. By analogy 
with (11), one now tries to choose a set of gains so that 


c = K(y) 


A r-A -i 


X 


— 5 


A 

For a given set of control surface locations, K(y;) is chosen to mini- 
mize some weighted, quadratic measure of the responses (e.g. , their 
r.m.s. values) over all time. 

But if y is fixed, the optimum is constrained in comparison to 
the values it might attain for some other y. The surface location prob- 
lem is to find the value or values of y which yield the best performance 
of the vehicle. The modified Newton -Raph son search procedure de- 
rived in the following section is a simple algorithm designed to ac- 
complish this complicated task. 



SECTION II 


THE MATHEMATICAL PROBLEM 

This segment of the thesis solves the mathematical version of 
the control surface location problem motivated by the discussion in 
Section I. A rigorous mathematical formulation of the problem is 
stated in the first subsection; the remaining subsections deal with the 
estimation, control, and search problems, respectively. Detailed 
calculations have been relegated to appendices. 

A. PROBLEM STATEMENT 

The linear, time -invariant system to be investigated in this 
thesis takes the following form: 

t) = F(y) + Q£.(x.< t) + Cri(t) (1) 

where x Ly> t) is the state of the (distributed) system 

.c (y,t) are control signals to the aerodynamic surfaces 
jq (t) is a vector of white noise system distrubances . 
y is a vector of control surface locations 


and x, jrj, and y are of dimension (n Q + n^), n c> 
spectively. 


n and n , re- 
T] y’ 


The vector t](t) is zero -mean with covariance matrix 


F(x(t)r{(T)) = N6(t-T) 

where N is a positive definite and symmetric matrix. Because N is 
time -invariant and the system itself is time-invariant, all time - 
varying signals will become very nearly ergodic shortly after t=0, re- 
gardless of initial conditions (which are left unspecified for this reason) 
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and the problem will become an algebraic one in the covariance 
matrices of the state variables. 

The state x (;y, t) of the system and the noise (t) are assumed 
to be of the form 


X(x,t) = 


; n(t) = 

"a o ro" 




_a 2 ( t )_ 


where x ,x,, n , -n-> are of dimension n , n, , n , and n , re - 
— o — 1 — 1 o -*£ o 1 o m 

spectively, and the plant and noise matrices are decomposed as follows: 


FCZ) = 

-F 0 G o (y)- 

; G = 

" 0" 

; c = 

1 

lo 

0 

lo 

1 

; N = 

'2 

0 

lo 
1 


-4 *1 - 


l 


_0 0 


^ ~ 2 _ 


The matrix G (y) is assumed to be continuous and differentiable in 
— o — 

each dimension of jy. (I is the identity matrix . ) 

The measurements of the state take the form 

= Ax(j,t)+ Bi](t) {2) 

In accordance with the above partitioning, A and B are written: 

A = [ A q 0] ; B = [0 B q ] 

B is assumed to be of full rank, n < n , for cases of interest, 
o m o 

The responses to be minimized contain state and control vari- 
ables: 

r(j,t) = H(i)x(x,t) + P(.y;)c(x,t) (3) 

The dimension of r (y , t) is and H(y) is partitioned as follows: 

H(y) = [Ho^yJFj 

Again, D(y;) is assumed continuous and differentiable in yr. 

The objective is to find a linear time -invariant control law of 


the form 
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c(£,t) =K(J)x(x,t) 


x(Y,t) = 


A , 

x 2 (Z> t) 


(4) 


where x (v , t) is the conditional minimum- variance unbiased linear 
1 o 

A 

estimate of given m(y, t), and K(^r) is chosen so as to mini - 

• mize the following quadratic criterion in the responses 


d(j) 


lim 

T-~oo 


T 

\ f t))dt 

0 


(5) 


where Q is a positive definite symmetric weighting matrix on the re- 
sponses. The surface locations, y, are chosen so that 

2* (y) < J*(y) for ally; in the domain 
where y is defined (assumed to be convex and s imply - c onne cted) , 
and where J*(y) is the minimum value' of the cost functional for a 
fixed y. 

In summary, the solution of the problem consists of three steps: 


(1) 

Find 


(2) 

Find 

* 

K(y) for a specified y 

(3) 

Find 

A 

2L 


The following subsections treat these problems in order as the esti- 
mation, control, and search problems. For the sake of legibility all 
y and t dependences will be assumed implicitly in the following 
derivations . 


B. SOLUTION OF THE ESTIMATION PROBLEM 

Kalman^ has shown the existence of a linear minimum- variance 

unbiased estimator x of the state x. The estimation problem may be 

^ Kalman, R, E., and Bucy, R. S., "New Results in Linear Filtering 
and Prediction Theory, 11 J ournal of Basic E ngine ering , March, 

1961, pp. 95-108. 
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described as follows: 

Given the system (1) and the measurements (2),, find the esti- 
mator of the form 

t 

X(X> t) = f h(X> t,T) m(y, T)dT (6) 

0 

which minimizes the criterion 

J(x) = TR [ E{(x-x)'(x-x) }] (7) 

assuming the control to be of the form £ = K(y)x(3r,t) . 

Make the definitions 


<i> 

X = X - 

A 

X 

= 

estimation error 

(ii) 

E{xx'} 


X 

- state second moment matrix 

(iii) 

E{x£'} 

= 

A 

X 

= state estimate second moment matrix 

(iv) 

E{x x 1 } 


X 

= error covariance matrix 


One then obtains the well-known equations for the Kalman -Bucy optimal 
filter L(t,T): 

9L(t,T)/8t = [F+GK-L(t,t)A]L(t,T) (8) 

where L(t,t) = [XA 1 + CNB'] (BNB')" 1 (9) 

and X is the solution of the steady-state error -covariance equation: 

X - (F - CNB'CBNB’r^X + ^(F-CNB’CBNB 1 ) -1 ^)' 

- XA * (.BNB 'f 1 AX + C(N-NB , (BNB , ) _1 BN)C ' = 0 (10) 

As a side result, when this filter is used, the estimation error has the 
property that 

E{x(y, t)*' (y.t)} • = 0 (11) 

which is a result of the orthogonal projection lemma.^ This in turn im- 
plies that 

A Kalman, R. E . , and Bucy, R .S . , "New Results in Linear Filtering and 
Prediction Theory, " Journal of Basic En ginee ring . Mar ch, 196 1, pp, 95 - 108 . 
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X = X+X ( PZ) 

Coming back to Eq. 10 and using the partitioned matrices given in 

» 

Section A, one sees that X must be of the form 



0 

0 


where X q satisfies the matrix Riccati equation: 


F X + X F' - X A' (B N^B^A X + ChUC' = 0 

— O ” O —O—O cr~ O O — 2— O -0—0 


(13) 


This equation is independent of y, a fact of considerable importance, 


it will turn out. In Appendix A it is shown that the state covariance 

a ~ 


equation can be expressed in terms of X(Y^) an d a constant 

matrix. The net effect of this is that, via (12) there is a separation of 
cont ro l and estimation problems and via (13) there is furthermore a 
separation of search and estimation problems? Equation 13 maybe 
solved once and for all and Eqs. 8 and 9 need only be evaluated after 
y and K(jy, t) are found--they play no further role in the problem so- 
lution. These side benefits of optimality show that use of the Kalman 
filter is justified not only by the fact that it serves as a standard of 
optimality but also by the great savings it yields in computational effort. 


C. SOLUTION OF THE CONTROL PROBLEM 

First, Eq . 5 is transformed in order to put the criterion in a 
form which is easier to manipulate. 

T 

J{j) = lim ~ f (_r r Q.r)dt 
T — oo 1 J 


0 


n n 
r r 


= V VQ.. lim 

A Af 1J T— oo 

1=1 j=l 


J. 

i/( v )dt 
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n n 
r r 


*2 s Q. . E {r^r.}, via the ergodicity assumption 

i=l j-1 


n n 
r r 


= ^jQ.j R.j, where R = is the response 

second moment matrix 


J(X> = tr[qr] 


(14) 


But R = E{_r r’} = E{(H x + D c)(x'H' + c’D')} 

= E{[ (H + DK)x iHx] [x^'H 1 + x'(H + DK) '] } 

= (H + DK)X(H + DK)' + HXH' 

Thus the criterion may be written 

J tx) = TR[ X (H + DK) 'Q (H+DK) + XH'QH] ( 1 5) 

/\ 

In Appendix B, the optimal control law K is derived via the matrix 

2 

minimum principle. The result is: 

K = -(P'QD)' 1 (D'QH + G'S) (16) 

where £> is the solution of the co state equation, 

(F - G( D ' QD)~ ^D ' QH) 1 S + S_( F -G(D 1 QD ) -1 D 'QH) -SG( D' QD ) 1 G 'S. 

+ H , (Q-QP(P , QP)' 1 P'Q)H = o (17) 

The state covariance equation for the optimally -controlled system is 


(F + GK)X + X (F + GK)' + FX + XF' + CN C 1 = 0 (18) 

A 

The Riccati equation (17) is solved first, yielding K which allows one 

A 

to solve (18) for X* All of this is carried out at each location y to be 
examined. The search problem dictates which y are selected. 


^Athans, Michael, n The Matrix Minimum Principle, 11 Information an d 
Co ntrol , Vol. 11, No. 5, pp. 592-605. 
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D. SOLUTION OF THE SEARCH PROBLEM 


Consider a s calar- valued function of a vector, J(y) which is 
everywhere continuous and twice -differentiable . If y^ * s a local 
relative minimum of J(y), the function may be expanded about this 
value to give 


J(x) = J(z Q ) + 

th 

where J^^) is a vector whose i element (i=l, . .. ,n^.) is 

hi^o) = [SJW/Syj]^ 
th 

And ^(y ) i- s ^matrix with i,j entry (i, j= 1 ... n^} : 




(19) 


( 20 ) 


( 21 ) 


This M element differentiation' 1 notation will be used throughout 
the thesis in a similar fashion for expansion of matrices about any 
nominal value of y. The subscripts (li) and (2ij) on matrices denote 
differentiation of each element by y^ and y^ and y^, respectively. 
For deterministic matrices (e.g., F(y), G 0 (y), etc.) the variations 

of X * are independent in each dimension, whereas for stochastic vari- 


ables (S(j), K(y), X(y» 
so that if y^ is varied, 
Returning to (19), 


the variations are taken mut atis mut andis 

y^, j/ i are adjusted to maintain optimal control. 

the conditions for a minimum v are: 

o 



(i) 

Ii(x 0 ), =o, 

(22) 

and 

(ii) 

(jz-x 0 ) , I 2 (Y 0 )(y-.Y 0 )> o 

(23) 


for all y in the neighborhood of y^ (this implies that J^(x q ) must be 

k 

positive semidefinite) . Assume that one has selected a point y which 

he presumes to be in the neighborhood of y^ and he wants to determine 

k k k k 

a Ay such that y 4 - Ay = y . One method of determining Ay is 



-21 - 


the Newton-Raphs on method, which proceeds as follows: 

Expand .Jj(x o ) about y Q : 

L^Z) ~ 1^) ~ J 2 (X 0 }Ay k 

~ IiCZq) - I 2 (x k ) A x k 

k k 

since _J-> (y ) differs from by second order in Ay . But 

_J-^{y o ) =0^ so one can write 

Ay k r- [J 2 (x k )] _1 ( 24 ) 

k 

providing the inverse exists, which will occur if y is in a neighbor - 

k 

hood sufficiently close to y^. Positive definiteness of J^(y ) can be 
shown a sufficient condition for convergence of the algorithm 

k+1 k , , k 

X = X + A X ( 25 ) 

to the desired value of y . The algorithm will give rapid err or- free 
convergence to local minima provided that is smooth enough 

near y^ to assure a large neighborhood of convergence. If convergence 

conditions are not satisfied one must use a gradient or other method to 
k 

bring y into the relevant neighborhood, resulting in a so-called 
"modified 11 Newton-Raphs on search. In aircraft problems, scrutiny 
of Section III will reveal that the requisite smoothness properties may 
be directly related to the smoothness of the mode slope derivates, 
which are in turn related to the smoothness of the elastic modulus as 
a function of position on the aircraft. Further discussion of this point 
is postponed to Section V, 

The search procedure used in this thesis is constrained in that 

k 

one demands an optimal controller be used at each y ; this reduces the 
criterion J from a functional in the two quantities y and K to a func- 

A 

tional in y alone by setting K= K(y) . The constraint equations are 
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derived in Appendix C, while the criterion variations are found in 
Appendix D; the results are: 

J = -TR[S 1 -(FX + XF 1 + CNC *)] (26) 

and J 2 . . = 2TR{[ (SF 2 . ,+S^j. + S^.+IH+DK) 'S(H 2ij +5 2i j® 

+<B lj +£ 1 jS) l Q(H li +P ] .K) -K j'. (D 'QDlKj.] X}, (27) 

where is the solution of the equation: 

(F+GK) 'S u + S li (F+GK)+F 1 '. S+SF li +(H li +P 1 .K) 'Q(H+DK) 

+(H+Dig 'QCH^+D^K) = 0 (2 8) 

and K 1 . = - (D , QD)~ 1 f D^. Q(H+ DK )+D l Q(H li +D 1 .K) + G'S^] (29) 

A z k , 

The matrices S^, K and X. are the solutions of Eqs. 16-18 for x=y 
In summary, the complete solution of the problem is given by 
Eqs . 13, 16-18, and 26 -29 and the search algorithm defined by (2 4) 
and (2 5). The implementation of this procedure on the computer is 


discussed in Section III-B, 



SECTION III 


implementation of the solution technique 

The practicality of the proposed solution technique rests on two 
factors. First, it must be possible to describe the two matrix func- 
tions of a vector G (;y) and D(y) for a general aircraft in such a way 
that readily available vehicle data can be eas ily utilized in the search 
algorithm; otherwise each application would require too much 
"accounting 11 engineering to be worthwhile. Secondly, an efficient, 
fast, generally applicable computational scheme must be written in 
order that answers are economically obtainable. Both phases of im- 
plementation are described in this section. 

A. VEHICLE DESCRIPTION 

In modern applications, it is reasonable to assume the availa- 
bility of linearized equations of motion and linear modal analysis of 
flexure properties of large aircraft. The following discussion will 
demonstrate that this information, plus some evident parameters of the 
vehicle, are sufficient to determine all of the variables needed in the 
proposed optimization scheme. 

A. 1 Specification of Control Surface Types and Positions 

The aircraft is broken into four segments: (1) wing, (2) fuselage, 
(3) vertical tail, and (4) horizontal tail --this numbering sequence is 
followed throughout the thesis. AH coordinates are referred to the 
orthogonal axis system below. 
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(b) SIDE VIEW 

Figure 1 



The true aircraft axes need, not correspond to these axes, of course, 
but positions on the true axis system are defined by their projections 
onto these axes rather than by measurement of lengths directly along 
the true axes. Customary units are inches; the end point of an axis 
need not correspond to a zero co-ordinate value (for example, the zero 
of the fuselage axis may correspond to the tip of a nose probe rather 
than to the actual nose of the airplane) - -for this reason two vectors, 
YO= (y 10 ,y 20 , y 3Q , y 40 )' and YEND = (y le , y 2e >y 3e > y 4e >' define the 
dimensions of the aircraft. 



The vector y actually consists of the positions of each control 
surface which is to be optimally located. The following control sur- 
faces are permitted (again, numbering corresponds to the code used in 
the optimization program) :* 

(1) Combined ailerons 

(2) Differential ailerons 

(3) Combined leading edge slats 

(4) Differential leading edge slats 

(5) Combined spoilers 

(6) Differential spoilers 

(7) Combined elevators 

(8) Differential elevators 

( 9) Rudder 

Positive deflections should be defined in the positive axis directions 
according to rigid body coordinates for the given vehicle, e. g. , if Y 
were out the right wing and Z were down, differential aileron de- 
flection would be positive for right aileron down, etc. 

Each of the five types of control surfaces listed above is as- 
signed an axis which is assumed to be a straight line and is specified 
by two quantities, its angle of intersection with the fuselage (AXANG) 
in radians and the distance (AXX) of its root from the center of gravity 
{ forward being positive)'. The set of force producers to be varied and 
their positions at a given point in the computer program are specified 
by the following vectors and matrices: 

^ Hencefor.th,.-.u i nless otherwise noted, subscripts on the variable y 
refer to the surface types as denume rated here (where y^ is pre- 
sumed to be with reference to the proper body axis), and not to the 
body reference axis system (y^ , y-> , y^ , y^) outlined above. 



AXANG 


AXX 


(a) AILERONS, SPOILERS, OR LEADING EDGE SLATS (TOP VIEW ) 




(c) RUDDERS (SIDE VIEW ) 


F igu re 2 
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NFP(I), 1=1, NT --number of force producers of type I to be 
positioned, where I corresponds to the numbering above , e.g., if 
NFP(2) = 3, three sets of (differential) ailerons are to be positioned. 

If there are no control surfaces of a given type, or if their positions 
are to remain fixed, NFP for that type must be zero. 

LOC(I) , 1=1, NT --axis or ¥ 4 ) to which type I 

control surface is to be referred. If rudder position is to be optimized, 

for instance, LOC(9) = 3, since rudder .position is referred to the 

vertical axis, y^, and the rudder is control surface type 9- LOC(I) 

th 

may be zero if the I surface position isn’t optimized or remains 
fixed. To facilitate added types of control surfaces (e.g., canards) 
the dimension of NFP and LOC is specified as NT rather than 9 
throughout the program. 

YFP(I, J), 1=1, NT; J=l, NMAX --position of control surface 
of type I, referred to the axis as specified by LOC(I) . Aileron, leading edge 
flap and spoiler positions are referred to wing coordinates (y^), ele- 
vators to horizontal tail coordinates (y^), and rudders to vertical tail 
coordinates (y^). 

AXANG(I), 1=1, NT --angle between fuselage and axis along which 
the influence of surface type I is felt (see Appendix G.3) 

AXX(I) , 1=1, NT --distance of root of axis of control surface 
type I fore of the C.G. (AXX for the elevator would be negative, for 
example). This parameter is often called in aerodynamic liter- 

ature. 

This completes the specification of aircraft geometry and control 
surface locations. The control and response matrices, G q and D, 
depend on these positions in four important ways: ( 1 ) through the flexure 
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physics, (2) through diminished lift effectiveness of outboard surfaces, 
(3) through increasing torque arms of outboard surfaces, and (4) 
through modal acceleration contributions to responses . The following 
sections derive these effects from the basic physics of the aircraft in 
flight. The casual reader is advised to skip to the summary following 
Section III. A. 4 at this point. 

A. 2 Airframe Flexure Physics 

This section outlines some of the significant aspects of airframe 
flexure physics. The discussion is based on Bisplinghoff ! s statement 
of the problem^ and on discussions with C. R. Stone of the Honey- 
well Systems and Research staff. 

Since discussion will be focused on the y-dependence of the G o 
and D matrices in the plant and response equations, viz., the de- 
pendence of forces and responses on outboard control surface location, 
fuselage dynamics will not be treated in detail. Although mode shapes 
include fuselage displacement, the analysis of wing, horizontal and 
vertical tail dynamics is the major issue of aircraft flexure physics. 
Each of these segments may be modelled as a tapered slender beam 
fixed at one end to the aircraft body. For such a structure, two types 
of motions may be distinguished (a) torsion about the axis of the beam 
and (b) bending about an axis perpendicular to the beam in the plane 
of the wing surface (see Fig. 3). Taking y as the co-ordinate out the 
axis of the beam, the dominant characteristics of each flexure mode 
are described by stating torsional and vertical deflections of the beam 
structure as a function of y. For slender slightly -swept or unswept 
wings the pure mechanics of vibration are well approximated by 
Hooke’s Law and St. Venant's theory of torsion. Only when the 



BENDING (z) 


TORSION (0) 


F igure 3 


torsional elastic axis passes through the center of gravity of each 
chordwise wing segment, however, do torsional and bending motions 
become decoupled. Figure 4 portrays the physical variables of the 
flexure equations: 


t (y,0 


iQ (y/0 


z (y,t) 


F igure 4 
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modes F z (y,t) and T(y,t) arise due to structural coupling and 
damping- -these effects are explained in Appendix I; for simplicity 
they are neglected in the following derivation. When F and t are 
set to zero, the solutions are separable functions of space and time: 


z (y»t) = e(y) T(t) . (3) 

0(Y,t) = <l>(y) T(t) (4) 

The flexure equations become: 

m(y)(e(y)-s(y)<t>(y)}T(t)+(EI(y)e ,, (y)) ,, T(t) = 0 (5) 

(Iy(y)<Ny) ■' s (y)m(y)£(y»T(t) -{Gj(y)4>'(y)) ’T(t) = 0 (6) 

'* 2 

T r / j \ i n-i /f\ i 1 * “1 _ . . i • 1 1 . . 


If T(t) +a) T(t) =0, then the spatial equations may be written: 
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For specified boundary conditions at either end of the wing, * these 

equations have solutions for an infinite sequence co = oo^, . . of 

natural frequencies . The solution corresponding to u)=oo^, consisting 

th 

of ?(y)=^(y) and t l ) (y)-«|>j(y)> defines the i natural mode shape of the 
wing. 

The distributed force problem is solved by assuming the flexure 

oo 

pattern to be a linear combination of natural modes: z(y,t) = 2 £.(y)r|.(t) 

oo i=l 1 1 

and 0{y,t)=S <j>.(y)r|. (t) where r|.(t) is the normal co-ordinate 
i=l 1 1 1 


See Appendix F, 
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specifying displacement of the i natural mode. Inserting these pro- 
posed solutions into Eqs . 1 and 2 and using the orthogonality condition 
that 

S. 

M j 6 i j=/ [m(y)^.(y)^(y)-m(y)s(y)[ <j> i (y)^{y)+<})j(y)^ i (y)] +I y (y ) cf,. Cy)<j>j (y)] dy 

0 

where t is the length of the wing, it is possible to find equations for 
the r).(t): 

i 

M.q ,{t)+M.o)? r).(t) =f (F (y, t)£ (y)+T(y, t)<j> ,(y))dy = X (7) 

J J J J J ' Z J J J 

0 

In practical applications the engineer is usually given the £,.(y) and 
4>.(y), or similar quantities^ The Eqs. 7 for the r^t) are then ap- 
pended to the plant equations in the model. The computer program 
models the direct forces due to control surfaces (flaps) on a wing as 
point forces. This means that the force from a flap centered at y^ 

will be approximately written F (y, t)=F (t)6(y-y ) and r(y,t)=T (t)6(y-y ). 

ZOO 00 

The above force term thus becomes: 

x. (t) iXy ) + t (t)<My ) - (8) 

j O J O O T j w o 

To a first approximation the forces F q and the torques T q are due to 
the change in lift resulting from a unit deflection of a flap of unit 
length --this may be found from tabulated section lift curves as a 
function of y Q . The forces are related to the control signal u(t) 
through actuator dynamics and Kussner lift buildup of aerodynamic 
forces, as described in Section I. The y-dependence of the mode 
forcing terms requires storage of vertical and torsional natural mode 


See Appendix F 
See Appendix I. 
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shapes . Since the variational procedure requires the second deriva- 
tive of these terms, an accurate spline fit procedure is used to gener- 
ate second derivatives, which are then Fourier decomposed for com- 
puter storage. Mode shapes and their derivatives are easily generated 
to an accuracy of better than one percent using only two hundred 
storage locations per mode. 

A. 3 Forces Depending on Control Surface Position 

The force matrix G (y) picks up its dependence on control sur- 
face locations in three ways. (1) The rigid body equations of motion in- 
clude the direct lift and moment effects of surface deflections- -some of 
the moment arms are functions of control surface location; (2) the 
extent to which a given flexure mode acceleration is excited depends on 
where the force producers are located; (3) rigid body and flexure 
modes are generally not decoupled- -vertical acceleration at the center 
of gravity, for instance, may include rigid body acceleration plus that 
due to each modal acceleration. 

(1) Rigid Body Effects. ' 1 ' The rigid body equations of motion are as- 
sumed to be stated in terms of a subset of the coordinates, 

(a) z - vertical acceleration 

(b) 6 - pitch acceleration 

(c) y - side acceleration 

(d) ip - yaw acceleration 

* f 

(e) (j) — roll acceleration 

(f) x - forward acceleration 


❖ 


The linearized equations of motion are derived from first principles 
in Appendix E . 
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All of these quantities are assumed to be evaluated at the center of 
gravity. For cases {such as the C-5A) where these are not the co- 
ordinates of the available model, it is usually easy to form the re- 
quisite linear combinations of these coordinates or to adjust the co- 
efficient matrices to account for changes of scale. 

(a) Vertical acceleration : The total vertical force due to a unit 

combined aileron, spoiler, leading edge flap or elevator de- 
flection is to first approximation due to the altered airflow 
pattern over that section of wing or tail. The effect of a 
deflection is to change the mean camber line, hence the ef- 
fective angle of attack, hence the lift of that section. The 
magnitude of the resulting vertical force depends on (of 
course) the amount of deflection, the depth of the hinge line 
into the wing, and the local chord of the wing, which is a 
function of y^. Thus for a .flap of unit length and given per- 
centage of local chord in depth, the total lift per unit angular 
deflection decreases as the chord becomes shorter, i.e., 
as the flap is moved outboard (see Fig. 5, curve a). 


Z 



Figure 5 



-34- 


To allow for chord variation and flap design this curve is 
specified for wing and tail as a Fourier series. For in- 
stance, if one wanted the same force per unit deflection all 
the way out the wing (see curve b above), the flap area would 
have to be increased for outboard surfaces. The shape of 
this curve thus embodies the flap design assumptions built 
into the program--curve (a) is for "uniform" flap size, while 
curve (b) is for uniform authority per unit deflection (but 
nonuniform flap size). Curve (b) is likely to be more useful 
for study purposes (a point force of constant strength being 
easiest to visualize), while (a) may be more suited to trial 
design work. The shape of this curve will influence the opti- 
mization results, however. It is assumed that the same 
curve applies to all control surfaces on the wing (exept for 
a constant scale factor). Letting 5^ be the deflection of the 
control surface (see Section III. A. 1), the y-dependent 
part of z may be written 


z = 

y 


S C T (y • ) 8. + M . A . (9) 

i=l ,3 , 5, 7 i 

C L .(y.) = CLC(i) • C ;L (y.) 


CLC(i) = scale factor for i** 1 surface 

C^(y) = section lift curve for wing or tail (see 

above); is one function for i= 1,3, 5 

and another for i=7 


M.A. 


= modal accelerations at the center of gravity 
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(b) Pitch acc e leration : Because of swept wing structure, as a 
control surface is moved out the wing or tail it is farther aft 
of the center of gravity, producing (along with the vertical 
force effect) a pitching moment which is just the product of 
the vertical force times the length of the lever arm, as shown 
in Fig. 6 . 



Figure 6 


Hence 


9 = 2 C T (y.)[ AXX(i).-y./TAN(AXANG(i))] 6. + M. A. 

y i=l, 3,5,7 1 1 1 


( 10 ) 


% 


Note that in calculating the lever arm AXX is taken to be 

th 

negative if the root of the center of pressure axis of the i 
control surface is aft of the center of gravity. The same 
diagram holds for the horizontal tail, AXX(7) being greater 
in magnitude. 
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(c) Lateral acceleration : Since no set of control surfaces can 

exert a balanced couple of side forces, there are no direct 
effects of either combined or differential control surfaces 
on lateral accelr ations . Since lateral flexure mode ac- 
celerations at the center of gravity are small, they have 
been neglected --few practical flexure models include lateral 
accelerations of the various modes, 

(d) Yaw acceleration : All differential control surfaces give 

rise to yawing moments due to differential drag effects of up 
vs. down deflections ; rudder deflections also contribute , The 
expressions for the moment arms are shown in Fig. 7. 


yAXANG 


AXX 


C.G. 



CENTER OF 
PRESSURE AXIS 


L.A. = [ y. 2 + ( AXX(i) - y. COT (AXANG(i))) 2 ] 1/2 


(a) WING CONTROL SURFACES 



L.A. = AXX (8) - y g COT (AXANG(8)) 

(b) RUDDERS 

Figure 7 
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Thus the direct effects of control surfaces on yaw ac- 
celeration are: 

rp= S C (y)[AXX(i) 2 +y 2 (l+l/TAN 2 (AXANG(i}) 
y i=2,4,6 1 1 

-2y.AXX(i)/TAN(AXANG(i))] 

+ C n ^[AXX( 9) + y 9 /TAN(AXANG(9))] {11}' 

where 

C^. (y) = GSTOCK(^ , y.) • C T (y), is the drag coefficient 
i l-l 


and C T is obtained from the vertical acceleration 
l-l 

equation. 


C N = GSTOCK(^,y 9 ) 

Note that the section lift curve (C T ) for each type of 

l-l 

control surface is used in approximating the differential 

drag (Cq ), the constant coefficient term being stored in 
i 

a reference matrix GSTOCK (explained in Appendix G.3). 
Since the sweep of the horizontal tail is negligible, the 
effect of differential elevator drag has been assumed in- 
dependent of y. 

(e) R oll accelerat i on: Rudder deflection and all differential 

control surface deflections produce rolling* moments . For 
rudder and wing control surfaces, the moment arms are 
merely the y-coordinate values for those surfaces. The 
diagram for the elevator roll moment arm is shown in 
Fig. 8. Hence roll acceleration is given by 


«|> = 
y 


s c p <y i )y i 6 i +c p (y 8 )[yvT +y 8 ] 1/26 8 +M - A - 

i=2,4, 6,9 i ^8 B Vi 8 8 


(12) 
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where 

C p = GSTOCK(‘i, y.) - (y.^ i=2,4, 6 

1 i-I 

C p ^ = GSTOCK(^,y 9 ) 



[(YEND-YO )\ +y g 2 ] 1/2 


Figure 8 


(f) Forward ac c eleration : The drag due to both combined and 
differential surface deflections depends on local chord and 
is hence a function of y^- For a given flight condition and 
control surface, the drag coefficient is assumed to be in 
constant proportion to the lift coefficient. Hence the y- 
dependent portion of the x acceleration equation is written 

9 

k* = S (L (y . ) 6 . (13) 

y i=1 D ; 1 i 

C D (y.) = GSTOCK(x,y.) .C L ( yi ) 

i i 

where 

GSTOCK(x,y.) = constant coefficient for x equation, 

th 

i control surface 

(2) The forcing functions for the mode acceleration equations depend 
on force producer locations. This dependence is described in Section 
III. A. 2, where the following theoretical results are derived: 
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where 


where 


+ co* ti.(t) = {F aero +X.)/M. (14) 

th 

r| . = j modal acceleration 

£ = structural damping (. 01 03)* 

th 

to j - j modal frequency 

^AERO = aerod y namic forces consisting of combinations 

■X. X 
-1" 

of other A/C state variables (independent of y) 


2 ^- = aerodynamic accelerations due to control surface 

j 

displacements 
9 C L. (y i > 

= ~m/m U j (y i ) + d (yi )*.(y i )]6 1 

th th 

§ .(y.) = vertical displacement of j mode at position of i 
force producer 

d(y^)= moment arm between center of pressure and elastic 
axis at l force producer position (formerly re- 


ferred to as s(y)) 

<j>j(y.)=torsional displacement of A mode at position of l 
force producer 


M. = modal mass of j mode 
J 


Observe that structural damping (Z^to^rj^t)) and aerodynamic force 

terms have been added to the equations given in Section III. A. 2. In 
addition, the terms F q and T q have been replaced by (y.)6.(t) 

and C T (y.)d(y.)5.(t) respectively. 

_L_J - 1 11 

1 

The generalized masses are input as the constants ZZM(j). When the 


aircraft is approximated by a lumped parameter model with masses 
m(y^) = and lever arms s(y^.) = s^., I(y) becomes I(y^)=s^mj . 

The expression for becomes 

M. = Z2M(j) =Sm k (^(y k )-2s^ j (y k )? j (y k ) + 4*j (y k » 
k 

(Continued on bottom of next page) 
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= J [ m (y)^j (y) -2m(y)s(y)(^(y)^(y)) + I y (y)] dy 


aircraft 

M = mass of aircraft 
The y-de pendent forcing terms, X^/M^, are 
routine GD. The equation and geometry for d(y^) 


computed by a sub- 
are shown in Fig. 9. 



ELASTIC AXIS 


CENTER OF PRESSURE LINE 
FOR ith SURFACE 


a ] = AXANG(i) 
a 2 = EAXANG(i) 


d(y.) = EAXX(i) - AXX(i) - y. COT fy) - COT fy) 

‘ = EAXX(i) - AXX(r) - y! SIN (a^ - a 2 )/SIN (a^ SIN (c^) 


Figure 9 


Since the lift coefficients are computed for accelerations, it is 

necessary to multiply througli by the aircraft mass (M) to achieve 
dimensional agreement or equivalently, to divide M. by M. Thus 

2>J 

the constants Z2M(j) = M^/M, with units of (length) are the actual 

input parameters of the program. 

In the C-5A case, the mode shape data includes masses m^ (26 

for wing, 25 for fuselage, 11 for vertical tail, and 10 for horizontal 
tail) and modal deflections for eack y^ corresponding to for 


each mode . 
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In summary, the y -dependent terms for the mode acceleration equations 
are 

X. 


9 C lJ>V 


■’ij = Z2M(jj = S “SiTm ts j ( y i ) + d(y i >,t, j (y i )]6 i (15) 

y ^ ■t j 

(3) Mode Acceleration Cont r ibution to Rigid Body Variables : As 

noted in the rigid body equation descriptions (1), the z, cj> and 9 
equations contain contributions from the mode accelerations at the 
center of gravity.* The z equation simply picks up vertical acceler- 
ations { hq j) at that pointihence .the mode acceleration term of the 
G (jr) matrix is 

z (M.A.) = 2 £ (y )[X /Z2M(j)] (16) 

y j j j 

where ^CG = center gravity location and X^/ZZM(j) is given above. 

The 9 variable picks up the mode slopes along the fuselage (nonzero 
for symmetric modes only), hence 


J 2 


)[X./Z2M(j)] 


(17) 


y^ = fuselage axis 


• ♦ 

The roll ( (j) ) equation picks up mode slopes across the wings (nonzero 
for antisymmetric modes only), giving 

$y(M . A . ) = (y o ) X./Z2M(j) (18) 

J 

For reference and comment, the relevant equations are summarized 
here: 


See Appendix E. 
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( 1 ) Vertica l Acceleration : 

NI 

z = E [C T (y.)+ S gly )X/y )/Z2M(j)]6. (19) 

y i=l, 3, 5, 7 L i 1 j=l J ° G J 1 1 

(2) Pitch Acceleration: 


d = E [C T (y.)[AXX(i) -y./TAN(AXANG(i))] 

y -U* 1 1 

i=l,3, 5,7 1 

NI 8g. 

+ 2 0^ {y CG )X j (y i )/Z2M(j)] 6 i 
J-1 ^ 


(3) Side Acceleration : 

Y =0 (2D 

y 

(4) Yaw Acceleration : 

*i = S C T (y.)[AXX 2 (i) + yf{l+COT 2 {AXANG(i))) 

y i=2 ,4,6 D i L i-1 1 1 

-2y.AXX(i)COT(AXANG(i))3 1//2 6. 

+ C D [ AXX(9) + y 9 COT(AXANG(9))]6 9 (22) 


( 5) Roll Acceleration : 



x=2,4,6,9 


t c P c l. 

1 1-1 


NI 8 g.(y rr .) 

2 - 3 ^ --- g -X (y ) / z2M (j)]5 

4_i 8 yi j x 


+ [ c p C L ( y g) t T +y 8 

o ( 

NI 8|.(y r _) 

H- S 8 J „ Cq -X.(y 8 )/Z 2 M(yj] 6, 

j — X 1 

{ 6) Forward Acceleration: 


(23) 


* « 

x 

y 


. , S , , , C D C lA> 6 i + . ? ?, 8 C D C L. <^> 6 i 
1=1, 3,5,7 i l 1=2,4, 6,8 i l-l 


+ C 


D, 


( 24 ) 
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A . 4 Responses Depending on Control Surface Position 

The y-dependence of the response matrix, D{y), is due only to 
acceleration responses; stresses and stress rates depend only on 
mode shape values at given response stations and on the modal dis- 
placements, which are functions of time alone. Stresses, for in- 
stance, may be written: 


s(y ) = k • (bending moment at y ) 


(2 5) 


where 

k 

EI(y r ) 

^(t) 


s; '(y r > 


= k -EI(y r )S I ' '(y r ,t) 
i 

= k *EI(y r ) - E^ftjg.' «(Y r ) 
i 

= constant depending on stress station 

= stiffness modulus at response station 

tli 

= displacement of i mode 


th 

= mode slope derivative of i mode evaluated at response 

station y 
7 r 

Since the r|^(t) (and T|.(t), for stress rates) appear in the plant state 
vector (x j) and the are fixed, the stress equations have no y- 
dependent terms. 

Acceleration responses, however, maybe written 


2 (y r > = R.B. (y r )+se.(y r W (t) 


(26) 


where 


R.B. = rigid body terms 

i.k 

g.(y r ) = l n mode shape evaluated at response station y^ 

rj (t) = modal acceleration of i^ 1 mode 
^i 
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Since only vertical acceleration responses are to be considered, the 

rigid body variables Y, ip and x contribute nothing. If y^ is on the 

fuselage, c(> may also be ignored. For y^ on wings or tail 

» ♦ # ♦ 

R. B . (y ) = z (y ) + <j> y - 0 x 
y w r y 7 r Y y 7 r y r 

For y^ on the fuselage, 

R. B. (y ) = z (y ) - 6 x 

y w r' y w r y r 

where x is the distance of the response location fore of the center 
r 

of gravity (i.e., is negative for response locations aft of the center 

of gravity. For y^ on the fuselage this is = y^Q - Y x and for y^ 
on the other segments of the aircraft it is x^ = EAXX-y^- COT(EAXANG), 
where EAXX and EAXANG are for the appropriate segment. This as- 
sumes that response locations will be located on the elastic axis. 


The modal acceleration terms t} (t) are explained above . Note 

^i 

that once the y-dependent terms of G q have been computed, they may 
be used directly in forming the D matrix. If 2s is not measured in 
the same units as the q. and rigid body variables , the appropriate 
elements of DSTOCK can be used for conversion factors. 


A . 5 Summary of Section III. A 

The preceding derivations appear exceedingly complex .because 
they attempt to treat a large class of aircraft models and use standard 
aerodynamic conventions (which are intrinsically confusing to the non- 
initiated) so that the parameters of available linearized models may he 
conveniently inserted into the computer program. The approach and the 
physics, however, are basically simple. The approach is summarized 
as follows : the G q matrix has the force coefficients corresponding to 
a given equation of motion along each row and the forces generated by 
each type of control surface down each coaiumn (see Fig. 10). The 
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pre ceding derivation simply passes from element to element, giving 
formulas for the nonzero elements. 


Equations^ Control Surface 
of Motion''''-^ Type (c) 


° 

1 

2 

3 

4 

5 

6 

M 

8 

9 

z 

X 

0 

X 

□ 

B 

0 

B 

0 

0 

y 

0 

0 


3 

□ 

0 

0 

0 

0 

X 

X 

B 

□ 

□ 

□ 

□ 

B 

B 

B 

6 

X 

0 

X 



0 


0 

4 

f. 

0 

B 

0 

B 

0 

B 

0 

0 

3 

i 

0 

B 

0 

X 

0 

B 

0 

X 

B 

t i.(sym) 

X 

0 

X 

□ 

Q 

0 

□ 

0 

0 

JL 

Ti.(asym) 

0 

□ 

0 

□ 

0 

B 

0 

X 

B 


General G matrix* 

— o 

{x denotes ^-dependent term; 0 denotes no -dependence) 

Figure 10 

A similar procedure is followed for the D matrix. The physical 
effects are basically simple also. The y-dependent terms occur in 
(1) the forcing functions for the flexure modes, (2) the torque moment 
arms, and (3) accelerations (either at response locations or at the 
center of gravity) at any point on the aircraft. The first effect utilizes 
the mode shape data, while the second effect uses basic -aircraft 
geometry; the third effect combines both data inputs. In addition, de- 
pending on the way one wishes to specify flap design assumptions, 
there may be a y-dependent attenuation factor due to decreasing lift 
effectiveness of outboard control surfaces. Appendix G. 3 describes 
how the various terms are computed. 

See Section III.A.l for definition of control surface types. 
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B. COMPUTATIONAL SCHEME 

The computer program which implements the proposed method 
has been written for full-scale investigation of the location problem 
rather than demonstration purposes. The input data consists of plant 
parameters derived from the aircraft equations of motion, gust model 
parameters, Fourier series coefficients for mode slope derivatives, 
and certain parameters specifying aircraft dimensions and control 
surfaces to be studied. The model may include lateral and/or longi- 
tudinal equations of motion and symmetric and/or antisymmetric 
flexure modes. Provision is made for the following control surfaces: 
elevator, rudder, spoiler flaps, leading edge flaps and ailerons {all 
surfaces operated in combined and/or differential mode). Most types 
of state augmentation may be accomplished with ease--e. g. , actuator 
dynamics, unsteady aerodynamics, etc. 

The structure of the program consists of an executive module 
which forms the coefficients of the filter equation (Section II.B.l), 
the optimal control equations (Appendix B), and the variational 
equations (Appendices C and D), and a series of subprograms which 
perform evaluation and computational functions . The subroutines in- 
clude : 

(1) A data read -in subroutine (READIN) 

(2) A routine which solves the Riccati equation (POTTER) 

(3) A routine which solves the state covariance equation 
(STCOV) 

(4) Routines to form D(y;), G o (y) and their partial de- 
rivatives (GD, GPDP, and GPPDPP) 

A Fourier series evaluation routine (EVAL) 


(5) 
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(6) Eignevalue and msitrix inversion subroutines (ALLMAT, 
HSBG, ATEIG, TDINVR) 

B . 1 Main Program 

The block diagram of the executive module (shown in Fig. 11) 
bears out the natural separation of the problem given in Section II. B, 
viz., estimation, control, and search. The estimation and control 
problems are well-known and need no further comment here, as the 
implementation follows precisely the equations in II. B. At the search 
stage, the manipulations in the actual program are slightly more com- 
plicated than shown in Fig. 11 , since the matrix YFP(I, J) , whose 
indices denote type and number, respectively of a given control sur- 
face, must be rearranged as a vector in computing the variations 1^^. 

and J- . 

Z jk 

The search procedure used is a modified Newton-Raphson scheme. 
If the second variation matrix is positive definite, the Newton-Raphson 


algorithm is used, but if ^jk * S not P os ^ ve definite a gradient 
method is used. The sum of the absolute step sizes in each dimension 
is taken to be a constant (one which is just large enough to move 
rapidly into a region where J^jk P os ^ ve definite), but the ratios 
of the step sizes in each dimension are in the same ratio as the re- 
spective components of the gradient, . More sophisticated tech- 
niques may be used. 

Because of its generality the program as a whole requires a 
considerable amount of computer memory. The object program gener- 
ated by the Fortran IV source deck, including dimensioned arrays, 
requires about 25,000 words of storage for a 14th order example; ap- 
proximately 15, 000 of which are required for the program itself and 
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Main Program Flow Chart 



Figure 1 1 
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Main Program Flow Chart 


FORM K{y ) = -(D'QDf A (D T QH + G'Sfr )} 

_r ~ 1 ~ 

FORM: 

A = F + GK& 1 ) 

S = FX + XF' + CNC' 


SOLVE AX + XA' + S = 0 


J 


1 

OALL O 1 UvV 


FORM J(x) = TR[S(i)(FX + XF' + CNC’) + H’QHX] X (END OF CONTROL PROBLEM) 


IV AR = 0 


IVAR = IV AR + I 


FIND 3G/3y‘ , 8D/3y‘ 


CALL GPDP 


FORM: 


A = (F tGK(i ))' 

S = ((H + pK(2))'QD lk + S(i)G lk )K(i) + K&VtP^'QfH + DK(y'))+ G lk 'S{y 1 )) 

1 , , 1 

SOLVE AS lk + S lk A' + S = 0 CALL STCOV X 


FORM K lk .= -(D'Qg) 1 Q{ 1 Q(jftDKfar 1 )>tD l QD lfc g(y I ) + gik^ ( I )+ ^ lg ik ) 

1 , 

FIRST VARIATION COMPLETE’ 

FORM^J lk (x)=2TR[ iS(z)G lk +{ H+DK^ 1 ))^^)^^ 1 ^^ 1 )] | X 

I " 

II - 0 

1 

v 

n = n + l I 


ji = o 

1 — 

J L 

ji = JI + l 


Figure ll(Contu* ) 



















- 50 - 


Main Program Flow Chart 



Note: The variational equations given above refer to the C-5A trial 
model of Section IV (see IV -A -8), although the program is 
readily modified to handle the general case of Section II. 


Figure 1 l(Contd. ) 
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for arrays of fixed dimension. With the availability of third "gene ration 
computers, computation time is a more significant economic con- 
sideration than storage requirements, however. The program does 
attempt to optimize on computation time. This is accomplished through 
fast single -product multiplications, elimination of needless loops and 
precomputations of certain matrix products, and (primarily) by great 
efficiency in the major computational subroutines. As an example, 
on the IBM 3 60/40 computer three iterations of a fourteenth order 
case took 2.75 minutes of main processor time; such a computation 
involves solution of four Riccati and 12 state covariance differential 
equations, all of order fourteen, plus the criterion computations of the 
search algorithm. 

B . 2 Major Subrouti nes 

A significant portion of the time devoted to this thesis was used 
in writing and debugging computer programs; in addition, several 
standard subroutines were tried and compared on the practical C-5A 
example cited in Section IV, as well as on trial examples. As a by- 
product of this research, the author wishes to introduce three sub- 
routines which significantly improve computation times over many 
existing schemes for solving commonly-occurring control problems. 
The subroutines are written so that they can be readily incorporated 
into programs of other users. 

The first subroutine (called STCOV) computes solutions of state 
covariance differential equations with time -invariant parameters; 
the subroutine is written to compute steady-state solutions, but can 
be easily modified to give the solution at any given point in time,, if 
desired. The iterative algorithm was developed by G. Stein of the 
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Honeywell Systems and Research Center* and has been used extensively 
there, but is not widely known. The author 1 s subroutine, written pri- 
marily as an exercise, closely parallels those in use* at Honeywell, 
Appendix G. 3 contains a description of the algorithm and computer pro- 
gram. 

The second subroutine (called POTTER) implements a method 
first described by Potter** for solving the steady-state time “invariant 
Riccati equation. It is well -knownf that knowledge of the steady-state 
solution can be used to convert the differential equation to a state co - 
variance form, which may be solved for the full time -varying solution 
via the above program or another scheme. For the time -invariant 
case where steady-state solutions are required, the subroutine de- 
scribed in Appendix G.l has been demonstrated to be definitely more 
accurate and faster by a factor of approximately one system order 
than commonly available iterative methods. Such efficiencies result 
from the fact that Potter's method is algebraic rather than iterative. 

The third subroutine (called EVAL) is used to evaluate Fourier 
series for flexure mode shapes, but can be used in more general ap- 
plications. Given a set of n functions defined along r axes 
(f^ (y^. . . y ) , . . . ^(y^ . . . y^)) and the Fourier coefficients of the second 
derivatives of the functions along these axes the subroutine provides 
for evaluation of all functions, derivatives, or second derivatives at 

A still faster version of the algorithm has been developed for re- 
peated solutions; the original results are given in Ref . 16. 

See Reference 12. 

t 


See Reference 14. 
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any point specified by a co-ordinate value of one of the axes. The ad- 
vantages of the program, which is outlined in Appendix G.4, are very 
compact storage of mode shape, slope and slope derivative data, as 
well as rapid access to these quantities at any point on the aircraft. 

A similar algorithm is available for obtaining the requisite Fourier co- 
efficients, if these are not accessible to the engineer. 

The remaining major subroutines are less generally applicable. 
They specify the dependence of control forces and response's on con- 
trol surface position. These subroutines (GD, GPDP, GPPDPP) are 
documented in Appendix G.2. 



SECTION IV 


APPLICATION TO A TRIAL MODEL OF THE 
LOCKHEED C-5A TRANSPORT 

In order to demonstrate the feasibility of the proposed compu- 
tational scheme and to obtain insight into the significant factors of the 
surface placement problem, a mathematical trial model of the Lock- 
heed C-5A transport aircraft was developed. The model illustrates 
the essential aspects of the problem by inclusion of three signifi- 
cant flexure modes and only two ailerons (per wing) for which optimal 
locations are to be determined. 

The reader is cautioned regarding any interpretation he may make 
of the numerical data presented herein; the figures are expected to 
bear at most a casual resemblance to those which might be obtained in 
flight tests. First of all, the data used in developing the model was 
compiled during design of the aircraft --before the transport had ever 
been flown. It is well known among aeronautical engineers that pre- 
liminary mathematical model data such as this may differ quite signifi- 
cantly from models based on flight test data, especially with regard to 
flexure equations. Secondly, the trial model uses only a fraction of the 
full aircraft model --it does not correct for the effect of coupling to less 
significant structural modes, nor does it include lags for aerodynamic 
lift buildup due to control surface deflections, or for actuator dy- 
namics. Finally, the author readily acknowledges that the trial model 
is a one-man effort and no doubt suffers certain deficiencies from the 
fact that he is not an aerodynami cist and cculd not have been perfectly in- 
fallible throughout the extensive hand calculations required to compile 
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the final trial model. Qualitative tests have been applied to verify the 
basic behavior characteristics of the model (e.g., sign tests, etc.), 
and the roots for the free -aircraft model do agree quite closely with 
those supplied in the original data base report. The author alone ac- 
cepts full responsibility for the validity of the procedures described in 
this section and emphasizes that neither Lockheed nor Honeywell would 
use a model such as this for even trial design studies. 

A. THE SIMPLIFIED VEHICLE MODEL 

The essential features of this model follow closely those de- 

"L' 

scribed in the Data Base Repo rt submitted by Honeywell, Inc. to the 
Lockheed -Georgia Company in its Load Alleviation and Mode Stabili- 
zation (LAMS) study of the C-5A. The report uses basic aircraft 
model data supplied by Lockheed in a simulation study of the signifi- 
cant structural modes and suggests lower -order models which accurately 
depict the behavior of the vehicle for a given flight condition. Beyond 
this information, the surface location program uses mode shape data 
computed by Lockheed (Case 110500, prepared on 4/7/67) for vertical 
and torsional deflections of the first, third, and sixth symmetric 
structural modes for the given flight condition. 

A. 1 Vehicle Geometry and Fl ight Configuratio n 

In the notation of Section III.A.l, Figure l,the vehicle geometry 
is specified by the values of and y^^, plus the sweep angles and 

intersection points of wing and tail (AXANG and AXX) center of pressure 
axes and (EAXANG and EAXX) elastic axes. These values are as fol- 
lows: 

= [ 60, 277.5, 405.6, 20.3] (inches) 

ind = t 1300 ’ 2841 * 3 > 800.7, 385.5] 

Reference 6. 
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The center of gravity 


1369.315 inches on the 


axis and re- 


lative to this 


AXX(l) = 220.55 (inches) 

AXX(4) = -1469-51 
EAXX( 1 ) = 166.81 
EAXX(4)= -1464.0 
The angles are 

AXANG(l) = 1.157 (radians) 

AXANG(4) = 1.163 

EAXANG(l) = 1.230 

EAXANG(4)= 1.163 

The variables for the body axis (vertical tail) have been omitted, 

since no rudder control is used in the trial model. In short, the wing- 
span and length are about 230 feet, with the wings at a sweep of about 
25°, making the C-5A the largest aircraft aloft (see Fig. 1 of the 
Introduction). The true center of pressure and elastic axes actually 
bend back slightly about half-way out the wings, but they have been 
approximated as straight lines. 

The flight condition is a low -altitude' cruise with the following 
specifi cations: 

Gross Weight: 5 93,154 lbs. 

True Airspeed: 57 7 ft. /sec. (M=.533) 

z 

Dynamic Pressure: 292 lbs, /ft. 

Altitude: 10,000 ft. 

The moments of inertia (see Appendix E) are 

I n = 30.3 x 10 6 (slug-ft. 2 ) 

I 22 = 28.8 x 10 6 
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I 33 = 56.3 x 10 6 

I 13 = 2. 1 x 10 6 

All flaps are in an undeflected position. 

A. 2 Rigid Body Equations 

For the given flight condition the equations of motion decouple 
into longitudinal (x, z, 6) and lateral (y, <j>) triplets. By defining 

the angle of attack a~z-V^9 where V" o =577 ft. /sec. and the pitch 
rate 6 and noting that all vertical aerodynamic forces and torques 
about the y axis maybe expressed as functions of a and 9 , one 
finds that Eqs. 13 and 1 5 of Appendix E further decouple from the 
forward acceleration equation, giving: 

a = f z (a-)/m (1) 

* *T- 

(&)= tya, 9)/l zz (2) 

Thus if one is only interested in pitch and vertical acceleration he 
need only be concerned with two of the six equations of motion due to 
the decoupling of forces and accelerations. Since the most significant 
responses of the vehicle are vertical bending and torsion stresses of 
the wings, and vertical accelerations in the fuselage due to vertical 
gusts, it is sufficient to treat just these two equations in the trial 
model. Because these equations are coupled to the flexure mode ac- 
celeration equations, further manipulation of the rigid body equations 
(see Section I, Eq . 1) is postponed to the next section. 


A normalization constant <f>^ was used to convert 6 to units of 
inches so that numerical coefficient values are of uniform size. 
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A. 3 Flexure Mode Equations 

The original vehicle data includes 30 flexure modes, 15 of which 
are symmetric in y-> and 15 of which are antisymmetric in y^ . 
Obviously, only the symmetric modes can give rise to vertical and 
pitch accelerations, so only 15 are relevant in light of the preceding 
section. For the trial model to be valuable in terms of insight, it 
was desired to cut the number of modes to a minimum consistent with 
a meaningful flexure control problem, namely three. Fortunately the 
JLrAMS study had identified the first, third, and sixth modes as the most 
significant contributors to RMS stresses and accelerations. Although 
this predominance was not overwhelming and although all modes ex- 
hibited fairly strong intercoupling these modes were seclected. The 
first mode is a pure wing bending mode, while modes three and six 
show coupled wing and fuselage bending, with fore and aft fuselage 
positive for positive wing bending in mode three and negative for 
positive wing bending in mode six. The normal co-ordinates for 
these modes will be denoted r^, , and (corresponding to the 

r|j(t) of Section III. A. 3 (2)) to avoid confusion with the gust disturb- 
ances . 

In correspondence to the system model of Section I.B, define 


*i = 


0/cj, 


-2 


3 

+ 6 


3 . 
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The LAMS report shows that the equations on motion may be written 
in the form 



where n = index corresponding to lifting surfaces 

i = index corresponding to control surfaces 

til 

= Kussner lift growth function for the n surface 

til 

= Wagner lift growth function for n surface 

Zg= vertical gust acceleration 

6. = control surface deflections 
1 

and (*) denotes convolution. The precise details of this model are 
not relevant to the thesis, except to note that ^ and correspond 

to unit -gain first order delays with time constants which represent 
the delay- time for build-up of aerodynamic forces on each lifting sur- 
face due to a change in attitude. These delays could be incorporated 
in the states and but they were approximated as zero in the 

trial model (this would never be done in practical design models) to 
simplify the system. The matrices may be roughly identified as 
follows : 


M = 

A 
— o 


K = 




F 

-gn 

F . - 
ci - 


mass -inertia matrix 

matrix corresponding to terms such as V q Q which 

arise due to use of stability axes 

structural restoring forces (elastic modes) 

aerodynamic damping forces 

aerodynamic linear restoring forces 

gust forces 


control forces 
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When the Kussner delays are ignored, the equations of motion may be 
written: 


M 


a + (A C + 2A )a+(K+SA 2n)a ~ 2F {j^*Z }+F 6 

n n n & s 


(4) 


n g' 

From this point it is a standard procedure to rearrange the equations 
into state -variable form in and Because the control surface 

aerodynamics and actuator dynamics (represented by the K^) have 
been ignored, however, the equations take a slightly different form 
than Eqs . 1 and Z of Section I. B: 

[c] (5) 


1 

I* 

h— * 

1 

= 




Y 

M + 

YY 

lA 


_±2 \ 

2_ 


A 

ay 


In this case corresponds to the (fixed) control surface locations 

specified in the available model and corresponds to { ib i . 


A • 4 Control Surfaces 

The original model used the two large outboard ailerons (see 
Fig. 1, Introduction), combined inboard spoiler flaps, and combined 
elevator flaps for control. Spoilers are drag -creating flaps on the 
upper wing surface which are usually used for landing but are of some 
value in flexure control. For the trial model, it was decided to keep 
the elevators and spoilers fixed in position and to determine optimal 
locations for two sets of ailerons (in combined mode, of course) . 

Thus there is a "significant" location problem, since two 
ailerons, essentially, must be used to control three wing-bending 
modes. Estimating the various parameters for the variable -position 
ailerons is a fairly difficult task which may be resolved by either 
reference to handbooks (e.g., the U. S. Air Force Stability and Con - 


trol Handbook ) or by asking the question "If an aileron were located 
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at the wing- tip, what would the computer predict its force -coefficients 
to be in terms of the input parameters? 11 and then using the available 
data to ,T fit l! parameter values. The parameters to be estimated are 
the appropriate elements of the GST OCK matrix and the section lift 
curve scale factor CLC(l) (see Section III. A. 3, paragraph 1(a)). The 
control vector c takes the form 

s 
a 

a 

e 

where 6 ,8,8 , and 8 are combined spoiler, inboard aileron, 

outboard aileron, and elevator flap deflections (units were converted 
from inches to degrees to make the coefficients of uniform magnitude). 

A . 5 Wind Gust Filter a n d Gust Aerodynamics 

The trial model uses the well-known Press-Meadows^ wind 
filter and the same Kussner lift -buildup functions that were used in 

'V-'Y' 

the LAMS optimization study of the C-5A, The gust filter has a 
transfer function 

a g (s)/-n(s) = (N ) (s+V/^ L)/(s+V o /L) 2 (6) 

where a = gust "angle of attack M 

r| = white noise (E{^(t)} = 0) 

Press, H. , Meadows, M.T., and Hadlock, I., "A Re evaluation of 
Data on Atmospheric Turbulence and Airplane Gust Loads for Appli- 
cation in Spectral Calculation," NACA Report 1272, 1956. 

Honeywell Interoffice Correspondence, April 19, 196'8, W. A. 
GlassertoM. A. Bender, "LAMS Optimization Study of the C-5A." 
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V q = forward velocity (577 ft. /sec.) 

L = gust correlation length (1000 ft.) 

s - complex frequency 

In state-vector form this is written as 

d (t) = (-2V /L) a ft) + p, + (^ V /L)n(t) 


g 


g 


P 5 (t) = -(V o /L) Z a g (t) + (V o /^L)-n(t) 

The aerodynamic lift -buildup effect of a gust striking a surface is the 

sum of two first-order delays with different time constants; the 

delays are represented by p^ and p^ for the wing and p^ and for 

the tail and are all driven by a (t). Hence the gust filter and aero- 

§ 

dynamics take the form: 


* 3 = J 3 £3 + C 3 tjj 


= ^#4 + —4—3 


(7) 

( 8 ) 


where 


s 3 = 


P5 

a 

g 


x . = 

—4 


n i 


The mean square gust intensity was taken to be 

= 50 (feet/sec)^ 


representing a gust field of moderate intensity. 


A. 6 Sensor Equations 

The sensor configuration chosen for the trial model does not 
correspond to that of the LAMS study but is motivated by the results 
of the optimization study, which indicated that the optimal controller 
required used as much "lead” as it could get in the filter networks. 
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In order to get this lead and provide some discrimination among the 
modes, the model uses wing tip accelerometers (averaged signal = 
m^}, an accelerometer on the nose (m-,), and a rate gyro at the center 
of gravity (m^) . Hence 

= a^z + ( x QQ~ x -^-rj-,)b^( @/§ 2 ) c i(^i(y-\^ r ’p)* T *l * ^3 ^ r WT^" T *3 


+ ^y (y WT* T 6 ) + ^2 

m 2 = a l Z * + ^ 2 ^ + c l^l^ x N^ r l + ^3 ^ X N ^3 + ^ 6 ^ X N^ T 6 ^ * ^3 


m 3 = (^/ ^ 2^4 
where 

a^, b^, = conversion factors (such that are in g T s) 

X WT ” distance from nose to wing tip along the fuselage 

= distance from nose to center of gravity along the 
fuselage 

= 0 (position of nose on the fuselage) 

= mode shapes of first, third and sixth modes 
evaluated at sensor locations. 
cj >2 = conversion factor (see rigid body equations) 

r l 2 ,T ' 3 ,T l 4 = sensor noise for wing accelerometer, nose ac- 
celerometer, and rate gyro, respectively. 

By means of the state equations (5) it is possible to express the sensor 
signals in terms of x.j_> X 4 J the controls c_, and , and r \^ . 

In keeping with the comment on Eq. 6 of Section I-B, it is assumed 
that the measured signals have been modified by subtracting out the 
known contributions of the control The variables r\^ 9 , and -q 

are obstensibly (white) sensor noise, but are usually magnified in 
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order to account for modelling uncertainty. For the trial model, 

and r|^ were assumed to have RMS values of .02, giving a signal/ 
noise ratio very roughly on the order of unity. One could argue the 
merits of increasing these values still further. So the sensor 
equations maybe summarized as: 


— - ^ + ^2-2 + - 4-4 + -2^2 


whe re 


” m f 

and ti 2 = 

"V 

m 2 


^3 

m 3 


_ T1 4_ 


(9) 


A. 7 Response Equations 

The locations of stress responses (s.) and acceleration re- 
sponses (Aj) are shown on Fig. 1 of the Introduction. There are 
two stress locations each on the wing and body and one on the hori- 
zontal tail root. The structural significance of these locations is 
readily understood from the diagram. There is one acceleration re- 
sponse at the pilot's station and two in the cargo bay; the fourth 
location is included to minimize flexing of the tail segment. These 
are the same responses used in the LAMS study. In the Data Base Re- 
port the bending moment equations are given as: 



+ SF, {K * 6. } 
. — k. l x 1 
i l 


( 10 ) 


Corrected to eliminate control signal contributions (see comments 
under the IT Control Problem” paragraph of Section I.B). 
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The variables are defined as in Eq. 3, and 


b = 







I 


] 


the numerical subscripts fallowing those on the stress locations in 
Fig. 1, with M v n denoting the vertical and T, c n the chordwise 
moment contributions . Using the approximations in Section IV. A. 2 
these equations maybe rewritten as 


b + (2N ln ) a+ (2N 2n )a+ F k x 4 +Z k 6 (11) 

n n n i 

The stresses at each point are linear combinations of the bending 

moments , thus 

s. = Kb , s. = [ s 1 s £ s 3 s 4 s 5 ]' (12) 

where K depends on structural parameters of the vehicle. The ac- 
celeration response equations are constructed in precisely the same 
way that accelerometer signals were constructed in the previous 
section, yielding 



as a function of x^, and Hence the response equations 

may be written as 

X = H 1 x 1 + H 2 x 2 + H 4 2£ 4 + D(x o )S (13) 

corresponding to Eq. 7, Section I.B (with F,-=0). Again, the param- 
eters of D correspond to the fixed surface locations of the LAMS 
model and the parameters of the variable -location model must be com- 
puted by a technique similar to that outlined in Section IV. A. 3. 


A. 8 Summary 

Combining Eqs . 5, 7, 8, 9 and 13 the complete system can be 


written as : 
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or 

* = Zx + G(i) u + Ctj (14) 


with rn = [ A 4 £] x + [ 0^ B^] ji = Ax + Bjj (15) 

and £ = [H x if 2 -4~^- + 5CZ)£ = Hx + P(y)u (16) 

The numerical values for the coefficient matrices are given in Ap- 
pendix H. From the above discussion, the covariance matrix for the 
white noise disturbances is written as: 



The weighting matrix, Q, for the criterion function gives a weight of 
unity to all responses but s^ (on the wing) and (near the center of 
gravity) which were weighted at ten. It was anticipated that, in view of 
the mode shape data, s^ would be especially likely to pick up stresses 
due to bending modes, while would pick up primarily modal ac- 

celerations. 

The mathematical model formulated in Section II, it should be 
noted, does not precisely correspond with this one, due to the neglect 
of actuator delays and control surface aerodynamics. In this model 
F and H are independent of _y, while G depends on ^L- This fact does 
not affect the form of the answers for the estimation and control 
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problems, but it does change the variational equations. For this 
system the variational equations (corresponding to Eqs. 26-29 of 
Section II. D) become: 

J u = -TR[S U (FX+ XF' + CNC’)] (15) 

J,.. = 2TR[ -K'. D'QDK 1 .X+ K'(g' S..+ g' .S + G' S-. (16) 

2ij -lj hr- ~ v -lj-li — 2ij— -li-lj ' ' 

+ (H+D K) + D[.Q P lj K)X ] 

(F + GK)^ li +J u (F+G|^ + K'(Gj i S+ pj.Q(H+pK))+{SG li +(H+pK) , Qp i .)K= 0 

(17) 

K u = -(D» QD )~ 1 (D 1 li QDK+D l QD 1 .K+D 1 .QH+G 11 S4-G , S li ) (18) 

These were the equations implemented in solving the surface location 
problem for the trial model. 

B . OPTIMIZATION RESULTS 

The general purposes of this section are two-fold: (1) to answer 

the question posed in the introduction, i.e. , how much performance 
improvement can be achieved by choosing surface locations optimally, 
and (2) to point out aspects of the results which will be useful in formu- 
lating the strategy of the next chapter. Of course all results are 
stated in terms of the behavior of the trial model of the C-5A de- 
scribed in the preceding section,: to the extent that this aircraft is 
typical, the conclusions can be generalized. Because flexure contri- 
butions to responses are very significant, one may presume that 
flexure control and hence surface location become more important as 
the number of significant elastic modes increases --i . e . , the more 
flexible the vehicle. Surface location could thus be considered a 
more important problem on very flexible aircraft such as the B-52 
and less significant (except from the standpoint of passenger comfort) 
on most of the smaller commercial craft. Although a detailed analysis 
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of the assertion is beyond the scope of this thesis, one may also pre- 
sume (to the extent that structure physics remain linear) that surface 
location becomes perhaps crucially important in extreme flight con- 
ditions, where flexure modes (if uncontrolled) may become critically 
unstable due to nonlinear aerodynamic effects. 

B,. 1 Significance of the Control Surface Placement Problem 

The first problem to be adressed is the evaluation of the signi- 
ficance of surface placement itself. In order to evaluate roughly the 
influence that location has on the performance criterion, the root 
mean square values of some of the state variables and the perform- 
ance index J(y) is given for some sample locations in Table 1. 

Figure 1 displays the positions graphically- -the twelve cases span the 
possible domain of values for y^ and y^ fairly well. 

The first case differs from the others in that it represents a 
nonoptimal control system for the given values of y^ and y^ . The 
suboptimal controller was generated as a result of a programming 
error but happened to be a useful case for comparison; Note that this 
controller does far better than any of the flexure controllers insofar 
as it controls vertical velocity very tightly; it also happens to con- 
trol the variables x^ through x^ as well as- -if not better than--the 
optimal controllers. The nonoptimal controller does much worse 
than the others in controlling the sixth mode (x^ and Xg) . The 
weighting matrix given in the previous section, however, has been 
chosen to weight heavily those responses due primarily to mode six 
(from the mode shape diagrams of Appendix I, one sees that mode six 
is most significant at stress location two and acceleration station two). 
Indeed the value of J(]£) for Case I is 50 times as great as for the 



Table 1 


Significance of Flexure Control 


CASE: 

LOCATIONS: 
y i 


I 


II 


'1 


473.3 400.0 

y 2 886.5 800.0 

STANDARD DEVIATION OF: 
x, 2.081 54.21 


1 


III 

IV 

V 

VI 

VII 

VIII 

DC 

X 

XI 

XII 

274. 1 

297. 1 

473.3 

440.0 

3 82. 8 

330. 4 

473 .3 

659.0 

777. 7 

855.4 

874. 1 

697.1 

886. 5 

953.3 

996.1 

1044. 0 

886.5 

1077.0 

1 163.0 

1235.0 

54.95 

52.24 

55.33 

55.47 

55. 99 

56.06 

49. 98 

51. 90 

53.02 

53.81 


x. 


x. 


X 


x, 


X, 


'8 


.1921 .2060 .1660 .2940 .1609 .1471 .1484 .1496 .2675 .2621 .2480 .2347 

1.227 1.681 1.706 1.557 1.660 1.604 1.657 1.643 3.309 2.479 1.774 1.621 

1.009 .2939 .3363 .2407 .3467 .4040 .4239 .4541 .4122 .2212 .2225 .2465 

1.067 1.392 1.396 1.316 1.356 1.285 1.362 1.337 .7868 .8781 1.051 1.114 

.9043 .2183 .2103 .2062 .2049 .1903 .2012 .1976 .0586 .1157 .1418 .1555 

9.543 .8079 .8353 .7530 .8586 .8764 .8912 .9033 .7547 .8993 .8971 .8961 

.4436 .04920 .05119 .04495 .05280 .03575 .05585 .05726 .05901 .06550 .06550 .06496 


CRITERION VALUE; 

J (y) 971.4 27.369 28.110 2 5.861 28.469 28.949 29.621 30.065 

DEFINITION OF STATES AND VARIABLES: 


N. A. N. A. N. A. 


N. A. 


y l 

y 2 


X 


1 


X- 


X, 


X 


X, 


X , 




x. 


8 

J(Z> 


position of inboard aileron (inches) 
position of outboard aileron (inches) 

vertical velocity with respect to stability axes (inches/sec, ) 

pitch rate (degrees/sec.) 

deflection of mode 1 (inches) 

velocity of mode 1 (inches/sec,) 

deflection of mode 3 (inches) 

velocity of mode 3 (inches/sec.) 

deflection of mode 6 (inches) 

velocity of mode 6 (inches/sec,) 

TR[S(e)(F X+ XF 1 + CNC')] 
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optimal controllers, demonstrating the value of optimal flexure con- 
trol per se 1 . Having demonstrated the merits of optimal control, it 
is worthwhile to ask how much optimal performance may be improved 
by paying attention to control surface locations. As a measure of 
the effect of surface locations on optimal performance, one can take 
the ratio of some variables from the worst case (VIII) to the best 
case (IV). The ratios are found in Table 2: 

Table 2 

Standard Deviation Ratios: 

Worst and Best Cases 

Variable x. x_ x x . x _ x . x_, J 

1 23 45 678 

Ratio 1.072 .505 1.055 1.886 1.015 .958 1.200 1.274 1.162 

Taken toge the r, variables Xg and the criterion J seem to indi- 

cate a potential improvement of about 20 percent in the performance 
of the optimal controller due to optimal positioning of the control sur- 
faces.* There is, of course, no guarantee that surface positioning 
will be more or less successful than this estimate with respect to 
any other criteron. On the basis of intuition, it seems likely that 
optimal surface location will yield a significant performance improve- 
ment in cases where the criterion emphasizes a number of modes not 
too much greater than the number of control surfaces. A comparison 
of cases I, IV and VIII, for example, indicates that if Q were chosen 
to weight x^ and more heavily (say by letting Q(6,6) = 10., 


Note that the "worst 11 locations of Case VIII do not appear at all 
unreasonable upon cursory surveillance of the aircraft. 
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weighting the pilot's acelleration) control surface positioning might 
be less effective than in the above instance. Pitch rate (x^) control, 
especially, seems to experience a trade -off with control of the sixth 
mode. The trade-off in pitch rate appears to involve the last few 
"ounces" of improvement for the original criterion, whereas the 
trade -off between vertical velocity and control of mode six seems to 
occur on a large scale (compare I and IV), except for the last degree 
of improvement, when the goals become complementary. > 

In the interest of practiced application, a detailed discussion of 
the philosophy of choosing quadratic weights would be very valuable 
at this point, but the subject is far too large and intricate to be 
treated in this document. Two general conclusions, however, may 
be drawn from the C-5A results. First, there exist some applications 
and some reasonable performance criteria for which optimal surface 
locations yield a significant performance improvement. Secondly, 
such performance criteria will tend to emphasize responses (especially 
stresses) which are closely related to a limited numbe r of flexure 
modes (i.e., a number not too much greater than the number of con- 
trol surfaces). This second point is utilized in Section V. 

B. 2 Analysis of Results for the C-5A Trial Model 

The analysis of this subsection is directed toward the identifi- 
cation of significant trends which emerge from the computer results 
rather than a precise description of the detailed performance of the 
optimally-controlled aircraft. Publication of detailed results serves 
no purpose because neither the vehicle model, the criterion, nor the 
performance accurately represents the true aircraft. Instead, there 
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is an attempt to evaluate the search procedure and to find indications 
of the "strategy 11 of the optimal controller. 

The Newton- Raphs on search procedure did not prove to be very 
effective, since the "neighborhood" of the minimum appeared to be 
too small for the positive semidefiniteness of the second criterion 
variation to be attained. Although some of the locations approached 
this condition quite closely, both roots of never actually entered 

the right half plane together, so a gradient method had to be used. 

The failure of to become positive definite could indicate (1) very 

sharply localized minima, (2) a very "bumpy" criterion with many 
local minima, (3) a saddle point, and/or (4) boundary extrema. A 
detailed analysis of the results seems to indicate (3) m one region and 
(4) in a second region (see Fig. 1) . If one visualizes as a 

surface over Fig. 1, cases II-IV tend to indicate a "valley" leading 
to a (inverted) saddle point, while cases IX -XII point to a boundary 
extremum. These results are not inconsistent, as both give local 
minima . 

Search results indeed demonstrate this trend, as depicted in 
Fig. 2. Contrast runs one and three. The first run starts at (400, 800) 
and the inboard aileron moves further toward the fuselage, causing a 
decrease in the criterion. The third run begins at (473,886) and both 
ailerons move outboard until the outermost one reaches the end of the 
wing. 

Although the two runs were made with slightly different computer 
programs (run three, from a Honeywell computer, uses longer 
Fourier series for the mode shapes than run one, from the IBM sys- 
tem at M. I.T. ), they accurately demonstrate the initialization 
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problem involved. The initialization points happened to lie near a 
local maximum of the criterion and hence the ailerons could n slide 
down 11 either side of the hill, which they did. Reference to Table 1, 
incidentally, indicates that the minimum of run 1 is probably lower 
than the boundary extremum of run 3 . * 

Thus the initialization problem can cause severe consequences. 
The only way to avoid this problem is to develop some effective 
strategy for making a good estimate of initial locations . Not only 
should such a strategy help to locate global minima, but is should 
also bring the program into a region of positive definite much 

more rapidly than a random search. 

If one recognizes that the objective of the optimal controller is 
primarily control of the sixth mode, the surface placement strategies 
for the trial model should be fairly evident from the mode shapes, 
which are also shown in Fig. 2. Evidently, the two options are (1) use 
the ailerons to exert a force couple at locations 200 and 800 or else 
(2) exert the couple of forces at about 800 and 1300. One might even 
be led to suspect the superiority of locations 200 and 800, since forces 
at these locations are less likely to excite the other modes of vi- 
bration than a force couple at the end of the wing. In general, how- 
ever, as the number of significant modes and the number of control 
surfaces increases, the optimal strategy will become less and less 
evident on intuitive grounds. The different frequencies of the various 
modes and the varying degrees to which the gusts excite each mode 
will confound any simple theory based on mode shapes alone. 

Due to a programming error, run 2 had a reversed sign on the 
gradient search and should be read backwards to indicate proper 
behavior of the program. Viewing this run as a lead-in to run 3, 
one sees the strange (but reasonable) search patterns which can 
develop. 
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If one were designing a suboptimal controller to n do the job 11 
within loose constraints, it might be possible to form some crude 
generalizations about the manner in which various modes might be 
controlled and to go on to invent a n pure n strategy based only on 
evident principles (e.g., location of surfaces near mode shape maxima 
in a way such as to achieve an ability to discriminate in control of 
various modes). When one is designing optima l controllers, however, 
the considerations can be quite different from those involved in a 
rough guess. Note, for instance, how the rigid body variables 
and radically changed roles between the 11 rough guess 11 and n final 

ounce” in the C-5A example. At the ” rough guess” stage behavior of 
was very loosely constrained, while control of pitch rate, x^ was 
more tightly constrained; at the ’’final ounce” stage, better control 
of Xj tended to clos ely coincide with reductions in J, while control 
of was in actual conflict with the reduction of J. Because flexure 
control itself is a rather ’’fine tuning” sort of problem (as opposed to 
rigid body control), the optimal strategy must depend heavily on the 
exact nature of the criterion, or at very least on the effect of the 
gains gene rated by the criterion. There is still hope, however, of 
devising a strategy based on one (or a very few) solution(s) of the 
optimal control problem for a given criterion—a strategy which 
avoids the costly and difficult search procedure. The remainder of 
the thesis is directed toward development of such a strategy. 

In order for a placement strategy to be applicable, consideration 
must be limited to a reasonable number of flexure modes, since the 
analysis must obviously have something to do with the mode shapes. 

A diagram of the migration of the roots for the optimally- controlled 
C-5A model indicates why this is likely (see Fig. 3). The optimal 
control system damped the free aircraft rigid body response 
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cons ider ably and also modes one and three. The important reduction 
of mode six responses, however, was accomplished almost as much 
by a frequency shift as by simple damping. At the higher frequency, 
less energy is put directly into the mode since both the gust cor- 
relation filter and the Kussner delays act as low-pass filters . As- 
suming that structural coupling of the modes is not dominant in the 
aircraft physics or aerodynamics, the effect of higher frequency 
modes cannot be strong, because a limited amount of energy is 
coupled into them by the wind gusts . 

As evidenced by the discussion above, most of the behavior of 
the search program for the C-5A trial model can be explained by 
direct reference to the flexure physics, so no note has been made of 
the actual values of the c rite rion variations . For reference purposes, 
the first variations for cases two through twelve are given in Table 3. 

Minimization of the first variation of the criterion (provided 
that second variation conditions are satisfied) is the true mathematical 
goal of the solution technique. In higher order problems, one is 
forced to rely increasingly on mode slopes rather than the mode 
shapes themselves, because mode slopes (upon which depends) 

give a better picture of the trade-offs between different locations. Be- 
cause it yields equations (i.e,, = J3) rather than the inequalities 

of global investigations, the strategy of driving to zero seems a 
good one to pursue. This approach also gives rise to a well-defined 
surface location procedure for the general case, as shown in Section V. 

One final trend in the results of the C-5A study is worthy of note. 
Once one of the two basic control strategies had been chosen, the 
optimal costate matrix and the optimal gains changed relatively 



Table 3 



AILERON LOCATIONS: 


y l 

400. 0 

274. 1 

i-H 

c- 

o 

CM 

473.3 

440. 0 

3 82.8 

330.4 

473.3 

659.0 

777.7 

855.4 

y 2 

800. 0 

874. 1 

697. 1 

886.5 

953.3 

996.1 

1044. 0 

886. 5 

1077.0 

1163.0 

1235. 

CRITERION VARIATIONS: 










9J/3 Yj 

.006458 

-.001432 

.09866 

. 002623 

.003841 

.005062 

. 006771 

. 04541 

-.00993 7 

-.02761 

-.03457 

9j/8y 2 

-.003 802 

. 01101 

-.1642 

-.005271 

- . 002873 

-.00460 

-.002615 

-. 1606 

-.06712 

-.03337 

-.01583 
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little as the search procedure progressed. The greatest variations 
appeared in the G (;jr) matrix and in the state covariance, although 
(see Table 1) even the state covariance did not change radically. 
These results of the C-5A example are perhaps the most important 
facts to note in developing a surface placement strategy. 



SECTION V 


A CONTROL SURFACE LOCATION STRATEGY 

The apparent complexity of the surface placement problem is the 
basic impediment to its solution. A wind gust strikes the aircraft, 
exciting both the rigid body and its elastic modes. The energy im- 
parted to the vehicle is dissipated through aerodynamic damping, 
structural damping, and transfer (via mode- coupling) to high frequency 
modes. The relationship between these forms of excitation and dissi- 
pation is very complex. The object of flexure control, generally, is 
to minimize the amount of energy dissipated in structural damping 
(since structural damping is ultimately a result of the heating and 
friction which occurs when the vehicle is deformed), especially at 
certain critical points on the aircraft, as this energy dissipation 
eventually gives rise to dangerous elastic fatigue. The optimal con- 
troller must disperse the gust energy in a way such as to avoid dissi- 
pation at these critical points. It can accomplish this (1) by "antici- 
pating 11 the gusts and keeping the aircraft in such a position that it 
abs orbs the least gust energy, (Z) by dissipating the energy in the 
viscous air flow (e.g., spoilers or rigid body motion) , or (3) by forcing 
the energy into elastic modes giving rise to noncritical responses. The 
C-5A flexure controller does all of these. The trade-offs involved in 
these forms of dissipation are very difficult to quantify. 

In order to solve the surface placement problem by a process in- 
volving less work (and more insight) than an iterative computer simu- 
lation of the aircraft, one must simplify the problem so that the 
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fundamentals are bared* The physics are simplified in Section V. A* 
Secondly, one must simplify the mathematics, since some ability to 
evaluate solutions globally is desirable, while it is impossible to 
carry the explicit y- dependence of all matrices through the entire 
problem. This is accomplished in Section V . B. Finally, it is neces- 
sary to display the alternatives in such a way that they are easy to 
evaluate, a task carried through in Section V.C. A demonstration of 
the strategy for the C-5A is given in the concluding subsection. 

A. SIMPLIFICATION OF STRUCTURE PHYSICS 

The following simplifications of the problem are introduced at 
this point 

(1) Assume point forces with magnitude independent of location, 

i.e., use section lift curves C T (y.) = CONSTANT. (See curve (a), 

l 

Fig. 5, of Section III. A. 3). 

(2) Linearize expressions for all rigid body moment arms. 

(3) Assume only stress (or stress rate) responses are to be 
used, so that the D matrix becomes independent of y, and Dj.=0, 
D-,.j=0 for all i,j = 1.. ,-ry. Actually any responses depending on the 
state x and not on x are permissible. 

(4) Neglect actuator dynamics and use the equations presented 
for the trial model (Section IV. A. 8, Eqs. 14-18). This assumption is 
not essential, but merely simplifies many of the variational equations. 

Now consider the effects embodied in the G(y) matrix. First, 
there are the forcing functions for the mode accelerations. Using 


Assumptions (3) and (4) are readily discarded, while assumption (1) 
may be relaxed by multiplying all mode shapes by C.^(y) suid adding 
some other minor modifications 



assumption (1) these maybe written (see Eq. 15, Section III. A. 3): 


/ 

‘*i = X C ij [ l i (y j ) + d(y j ) ' l> i (y i )3 6 j 

^ i=l 


( 1 ) 


where c.^ - ig independent of y.. In the spirit of the 

summary of Section III, a shorthand notation for the elements of the 
G matrix is now introduced. Each element is denoted by G(EQ; s), 
where EQ is a symbol describing the force equation and s is a sym- 
bol denoting the type of control surface. Using this notation, the co- 
efficients corresponding to (1) may be written as 

G(\, j) = c..(§.(y.) + d(y.)<j>.<y.)f (2) 

It is also advantageous (in light of assumption (3)) to define the modified- 
mode shapes . fh(y), as follows: 

^(y) = £.(y) + d(y)(j>.(y) (3) 

With the convention that 


P Q (y) = y 

Thus (2) becomes 

G(\,j) = c..(3.(y.) 

One can similarly write the rigid body forcing terms as 

NI 

G(M k ,j) = g k . + c fco .p o (y.) + £ c k ..p.(y.) 


i=l 


NI 


= *kj + Z °wj 


i=0 


(4) 


(5) 


( 6 ) 


.th 


Multiplication by 5., the i surface deflection, occurs when G(y) 
is multiplied by u, the vector of control settings. 
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where 


th 


x, = k rigid body acceleration (e . g. , z,X, 9) 


g 


, . = constant terms (e.g., C_ (y.) = C T in the 

kj Lt. j 

z equation) 


c. . . = coefficient of i** 1 mode acceleration due to i^ 

kij 

control surface in mode acceleration pick-up of 
th 


k rigid body equation 


.th 


c, . = moment arm coefficent of j control surface 
koj 

The evident course is now to define a matrix B(y), to be dis- 
tinguished from the constant JB matrix of the sensor equations, whose 
.. ..th 

(l, j) entry is 


where 


i = 0, . . ,,NI 
j = 1 , . . . , n 


(7) 


Finally, the G(y) matrix may be written as 


fi(2) = G o +G^Btz) 


(8) 


G q is a constant matrix containing terms such as g^. of Eq. 6 and 
also force coefficients for surfaces which are fixed in position. The 
matrix contains the coefficients c„ of (5) and c^.j of (6). 

With this notation it is also easy to represent the derivatives of 


G with respect to y^, G.^: 

dG(x) 

G.. = "5 

“lj 


= G 


8B(x) 


-1 9y- 


(9) 


th 

Note, however, that only the j column, b., of B(^;) depends on y., 

J J 

th 

so B, . has only the j column, 8b. /3y., nonzero. This fact motivates 
— lj J J 

the shorthand notation: 



-85- 


G^(X» = G x B^y) 


( 10 ) 


.th 


where the i column of Gj(y) contains the nonzero entries of G^., and 


.th 


the j column of B^(^r) contains the nonzero entries of B^, i.e 


9p. 

^i ( ^ )] ij = aF {y j ) 


(ii) 


9 Pi 


where = £j(Yj) + d , (y.)$ i (y,j) + d(yj)cj>!(yj) i/0 


9(3 

a-nd ^(y.) - 1 


U0 


is simply a matrix of the modified mode slopes evaluated at the 
various surface locations. 

Thus there exists s relatively simple expression for both the 
explicit y-dependence of G(y) and the ^-dependence of its first 
derivatives G ^ j (X ) # Furthermore, due to assumption (3), this com- 
prises the entire y-dependence of the problem. 


B. 


SIMPLIFICATION OF MATHEMATICS 


In the best possible situation, one would hope to be able to obtain 
an explicit expression for J(y) in terms of y itself, or in terms, say, 
of the mode shapes. It will be shown that this is not feasible. The 
next best situation would be to obtain the first variations of J(y) as 
explicit functions of the normal modes . Provided that' one is willing 
to make a few approximations, it will be shown that this is possible 
(although the optimal control problem must be solved at least once). 
Furthermore, the requisite expressions for JT^y), i = l. . . T] , are 
readily evaluated and put into a form where near -optimal locations 


become evident. 
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First consider the problem of obtaining an explicit expression 
for J(y). The most promising expression for J(;y) is given in Ap- 
pendix D: * 

J(X) = TR[_S{_£ ) (FX+ XF 1 + CNC') ] (12) 

This expression is promising because the matrix (FX+XF +CNC ) is 
independent of jy and depends only on the results of the state esti- 
mation problem, which can be solved independently of jy . The costate 
matrix p(;y), in addition, is the best one could hope for, in that K(;y) 

A 

and X(y) are derived from it, and would necessarily be more compli- 

A 

cated functions of jy (i.e., K(;y ) involves G (y)S(y) and X(jy) in- 
volves G(_y)K(;y)}. Unfortunately, however, S(_y) is a s olution of the 
Riccati equation 

(F -G(y) (D 'QP) _1 D 'QH) ‘S{y) + S(j) (F - G(y) (D , QP)“ 1 p 'QH) (13) 

-S(x)G(y)(P'Qpr 1 G , (Y)S(x)+H , (Q-Qp(p , Qp)" 1 p , Q)H = 0 

The third term of this equation involves, effectively, the product 

* which is a thorough mixture of mode shapes, moment arms, 
constants, and every combination of these three multiplied together . 
There is virtually no hope of carrying any n bookkeeping M procedure on 
the mode shapes through the solution of the Riccati equation . The linear 
equations which give hounds on the solution of the Riccati equation might 
yield relatively simple expressions in the mode shapes, but the bound- 
ing solutions are known to be very poor. Thus there appears to be 
very little hope of finding the criterion as an explicit function of _y, 
other than mapping it tediously, point by point. 

The next alternative is to seek the variations of the criterion as 
functions of One might well ask why this should offer any more 

hope than the evaluation of the criterion itself, as the first variations 
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have been seen to involve at least the solution of the control problem 
at each location, which is just the task that led to difficulties above, 
due to the squared terms in G(x )G r (x) • Some justification for using 
the first variation of the criterion can be derived from this very fact. 
Differentiation of the product terms can yield only: 

(1) zero 

(2) mode shapes 

(3) mode slopes 

(4) products of mode shapes and mode slopes 

But if the mode slopes are sufficiently strong functions of x as op- 
posed to the mode shapes , one may well be able to justify using some 
nominal values in place of the mode shapes, and thinking of the first 
variations as ne arly valid in a relatively large neighborhood of these 
nominal values . This argument, admittedly, is heuristic; neverthe- 
less, it does motivate an avenue of approach. 

The most obvious expressions for the first variations are derived 
directly from (12) above, i.e., 

JjjCz) = TR[_S li (y)(FX+ XF 1 + CNC 1 )] , UL..T, y (14) 

The difficulty with these expressions, however, is that the matrices 
_S ^i (jy ) are solutions of the state covariance type of equation given in 
Section IV, Eq „ 17. Although the functions _S^(x) have a linear 
dependence on the mode slopes, as can be seen by examining the row 
operations used in an algebraic solution of the system of equations, a 
great deal of IT bookkeeping M is involved, and this motivates one to look 
for a more useful representation of _J^(x)- Such an expression actually 
exists : 

J(i) = TR[S(i)(FX+ XF' + CNC')] 

= -2TR(S&)(F+ G<2) K(X)) t(x)] 


( 15 ) 
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Using the equations of Appendices B and C, modified for the 
dependence in the trial model equations, one finds that (for 
it is possible to write 

J u (2) = 2TRL.S&) G u (£) K(2 )X(y)] (16) 

C. THE SURFACE LOCATION STRATEGY 

In the spirit of the preceding discussion, one expects that the 
approximation 

\(Z) = 2TR[S( io ) G h ( Y ) K( Zo ) X( Zq )] (17) 

will be nearly valid in a relatively large region around y; o * The opti- 
mal control results of Section IV . B indeed indicate what sort of a 
region this will be. The C-5A study indicated two distinct control 
strategies, one using a force couple near the fuselage and another 
using a couple near the wing tips . In the y^, y^ plane, these strategies 
yielded optimal locations which were widely separated, and control 
laws which differed significantly between but not within the region of 
each strategy. Thus in the general case, one is prompted to think in 
terms of "regions" where different basic strategies prevail, the 
boundaries of these regions being (essentially) determined by the 
character of the weighting matrix, Q. As the number of modes or 
the number of control surfaces becomes greater, and/or the weighting 
of the modes implicit in Q becomes more evenhanded, the number 
of "regions" is likely to increase. One anticipates that an expression 
such as (17) will be nearly valid within the "region 11 of and per- 

haps somewhat beyond, with luck. 

Within each region, however, J^(y) has become fairly easy to 
evaluate explicitly, since, according to (9) , contains si mple 
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linea r terms in the modified mode slopes. It is a very easy matter to 
write a computer program which does the bookkeeping involved in 
evaluating (17). The expressions for the J^(^y) will be of the form: 

NI 

J u te) = ^(yj) = Y, a ij (8p j (y)/8y i )ia ij = a ij te o’ <18) 

j=o 

Following the heuristic notation of (10), one may write 

J^) = 2TR[_S{^ o )G 1 B 1 ( i ) K (l 0 )X(y Q )]^ (19) 

In summary, the following search strategy promises consider- 
ably more insight and more rapid convergence on a globally optimal 
set of control surface locations than a Tl blind n gradient (or Newton- 
Raphson) search: 

(1) Choose a reasonable initial set of surface locations, y , 
based, if possible, on a best guess of what the optimal controller will 
attempt to accomplish in view of the weighting matrix Q. 

(2) Solve the optimal control problem for _S(J o ), K(X Q ) 

A 

X(y Q ), as given in Appendix B. 

(3) Evaluate J^.(y) according to (17) above. 

(4) Using these functions, determine if y Q appears to lie in 
a large region where can be nearly driven to zero; if so 

(a) Choose a new set of locations and proceed with a 

modified Newton-Raphson search (or call the answer); if not 

(b) Determine the boundaries of the region near and 

make an educated guess at y on the basis of 


A similar expression for J^Cz) 


can be derived if desired. 
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(5) If 4(a) is true, return to (2) and proceed with a more exact 
local search; if 4(b) is true, repeat (2) -(4) in order to identify a 
better "region." 

D. DEMONSTRATION OF THE SURFACE LOCATION STRATEGY 
The above steps have been carried out for the locations 
jy Q =[ 297.1, 697. l] 1 of Case IV, Section IV . B . For the G(y; ) matrix 

G<d,j) = S (-3.25e i (y CG )/Z2M(i))(3.(y ) 

i=l,3,6 J 

G(*<9/4> 2 J) = 2 {-3.25e|(y CG )A 2 *Z2M(i))p (y ) 

i=l ,3,6 J 

+ (.791 x io" 4 /4> 2 )yj 

= (-3.25/Z2M(i))p.(y.) , i=l,3,6 

Inserting the relevant constants and taking the partial derivatives with 
respect to y ^ and y 2 (j = l,2), one obtains G^Cy): 

G^a , j) = 11.153p'(y.)-.96l23^(y.)-1.691p^y.) 

G . (d , j) = +.1304 -I- 16. 085p 1 '(y.) + 6.2554p'(y.)-4.3894p' (y.) 
Gj(t , j) = -103. 2 7(3 ’<y.) 

G^t^J) = -3 8. 449|3^ (y.) 

G^J) = -6.9592p;(y.) 

These equations are readily cast into the form of Eq. 10. Using the 
computer results for _S(X Q ), K(y Q ) and X(y Q ) and evaluating the ex- 
pressions (17), one obtains 

(J l ) 1 7 2. 1065 pjty^h . 06489 p3(y x )+ .4337 (3^) + .02826 
(J^ = “4. 13 98 pj(y 2 ) + 1 . 5l2 0p3(y 2 ) + . 7 9705pg(y 2 ) - . 03849 

These functions are sketched in Fig, 1, 
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The actual locations for case IV are denoted by (y°, y^) ; the 
gradient method generated a new set of locations (yj , y 2 ) which are 
shown in Fig. 2. According to the proposed surface location strategy, 
the optimal locations should be (y^ , y£) of Fig. 3, although some com- 
ments on these locations are in order. The first fact of note is that 
the bias terms (due to rigid body effects) are too large to allow the 
condition JT^=0^ to be met for any . This possibly indicates a model- 
ling error in one of the force coefficients, a problem which would have 
otherwise been very difficult to identify. Note that if one were to as- 
sume for the minute that rigid body torque are in effects were negligible, 

'VI 

these biases would become zero, yielding optimal locations yj — 280 
or 480 and y^ - 1060, which appear quite reasonable. Secondly, 
note the indication of the "regions' 1 . Beyond y^ w 600, (J^ blows up, 

'V' 

indicating a likelihood that (Jj) has passed out of its region of validity 
(indeed, for the other strategy, the optimization results indicate 
y^> 600 as in case XII) . Similarly, one observes that in the region 

) 2 experiences a maximum. This 
clearly indicates that y 2 is near a barrier between regions and that 
the behavior of the criterion may change considerably as y 2 is 
varied. These notions indicate the sort of global insight one might 
anticipate from the proposed scheme. Evidently, human cleverness 
has a considerable role to play here. 


of validity (y 2 « y°) for y 2 , (J^ 



SECTION VI 


SUGGESTIONS FOR FUTURE RESEARCH 

Many pos sibilities for future research into both mathematics and 
physics have been indicated throughout this work. Promising compu- 
tational methods give some hope of making formerly intractable prob- 
lems mere exercises in linear algebra. 

First, some basic investigations for distributed systems are in 
order. Some general convergence proofs must be developed in order 
that consideration of a limited number of flexure modes has rigorous 
justification. The energy considerations mentioned in the introduction 
to Section V, plus a bound on the amount of coupling between flexure 
modes should (with difficulty) establish this result. 

Next, a formalism must be developed in order to enable one to 
optimize on number as well as location of control surfaces. In this 
thesis, the number of surfaces is tacitly assumed to be prespecified. 

The form of the system might be generalized to allow sensor 
locations as free variables, in addition to control surface locations. 

In this event, decoupling of the estimation problem from the other 
phases of solution can no longer be accomplished, although use of 
optimal estimators may simplify the variational problem as it did in 
Section II. 

More application- o-riented research is definitely needed in the 
area of search techniques. The proposed search strategy should be 
fully investigated, as it offers some hope of being useful in obtaining 
global solutions . The current data is too scanty to allow adequate 
evaluation of how generally applicable the proposed strategy is. 
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Finally, there is the problem of simplifying the physics. Be- 
cause a great number of physical systems can be put into the form of 
the elasticity equations for the aircraft, a vigorous attempt should be 
made, in the spirit of Section V. A to obtain explicit general solutions 
for the distributed -force control problem and sensor location problem. 
Tensor notation, plus use of the close analogy between equations in 
space and time appear to be essential ingredients in the solution of 
these problems. 


rr - 

Research along these lines has recently been carried out* See 
Greenberg, S*, Optimal Pointwise Feedback Control of Distributed 
Parameter Systems , Doctoral Dissertation, Massachusetts Institute 
of Technology, June, 1969. 



SECTION VII 


SUMMARY 

The basic physical relationships involved in control of a flexible 
aircraft disturbed by random wind gusts were used in formulating the 
surface location problem as one in optimal control of a distributed 
system, using a limited number of point-force controllers. The three 
phases of this problem- -estimation, control, and surface placement -- 
were then solved by means of the matrix minimum principle and the 
calculus of variations. The variational equations were greatly simpli- 
fied and the order of the problem considerably reduced by the use of 
optimal controllers at each stage in the search for optimal surface 
locations. These simplifications, plus advanced computational tech- 
niques made the general solution of the problem practically feasible. 

Aircraft physics had to be investigated in great detail in order 
to obtain general equations expressing the distributed nature of the 
system exactly, A computer program was written which stored these 
equations and used them in solving the surface location problem for 
a general aircraft. 

This program was tested on a fourteenth order model of the 
Lockheed C-5A transport aircraft. As a guide for future applications, 
the derivation of the model parameters was carried through explicitly. 
The results of the optimization study were then analyzed in an at- 
tempt to recognize and develop a Strategy 11 for control surface 
placement. 
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A simple practical strategy was developed for systems with 
stress and stress rate responses and slightly simplified physics. This 
strategy, which was also verified for the C-5A model, has the two 
major advantages of (1) insight into the tradeoffs involved in surface 
location, and (2) partial insight into global solutions. 

The major contributions of this thesis are (1) a computationally 
feasible general solution to the surface placement problem, (2) a 
practical application of optimal flexure control to a trial model of the 
C-5A transport, (3) implementation of efficient computational tech- 
niques for solving state covariance and Riccati equations, and (4) 
recognition of the nature of the optimal solutions and presentation of 
a search strategy which makes use of the basic aircraft flexure physics. 



APPENDIX A 


DETERMINISTIC FORMULATION OF THE PLACEMENT. PROBLEM 

The state and state estimate differential equations are: 

x = F x + G Kx + Ci} 
x = (F + GK)x + L(m - Ax) 

= (F+ GK- ^x + L Ax + L Brj 


JQ = E 2 ~ + C 2^ 


Defining = 

[x , x] write 



z ~ 

F GK 

_Z + 

c 


LH F+GK-A 


LB 

_ 


The associated covariance differential equation is: 
Z = F Z + ZFl+C,NC; , where 


■ 2 — --2 — 2 — —2 
Z = E{zz'} =| 


A* 

X X 

A A 

X X 


E{xx’} E-fxx 1 } 

E(xx'| E{xx'} 

using the orthogonal projection lemma of Kalman. Extracting the X 
equation, 

X = F X + GKX + X F 1 + ±(G^' + CNC' (1) 

= QT+ GK}X+ X(F+ GK) ' + CNC ' - GKX-X K'G 1 


A A 


( 2 ) 

= (F+ GK)X-t-X(F+ GK ) 1 + FX+XF' + CNC 1 (3) 

since X = Efxx 1 } = E{(x + x)(x' +3c')) = Efxx'HEfxx 1 } = £ + X 


Similarly, X = X+ X , so 

X = (F + GK)X+ X (F+ GK) '+CNB ' (BNB ' ) _1 AX 


(4) 


4 X(CNB 1 (BNB l )~ 1 A) l +XA , ( BN B , r 1 AX+ CNB 1 (BNB') i BNC' 


-1 


[ u 
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Of these four formulations of the problem, expression 3 is most con- 
venient. The last three terms of this equation have no y-dependence 

and thus disappear in the variational equations. Similarly, the cri- 

A ~ 

terion function maybe separated so that J(X) = Jj(X,'y)+ J,{X) Thus 
also disappears from the variational equations. 



APPENDIX B 


DERIVATION OF THE OPTIMAL CONTROL (FIXED LOCATIONS) 
VIA THE MATRIX MINIMUM PRINCIPLE 

Form the Hamiltonian for the problem using the LaGrange multiplier, • 
S 

J = TR[ QR+ SX] 

= TR { Q[ H+ DK)X (H+ DK) ' + HXH 1 ] 

+S[ (F+GK)X+ X(F+ GK) 1 + FX+XF' + CNC 1 ] } 

A 

Taking as fixed, choose K, S and X to minimize, J: 

9j/9X = <H+ DK )‘Q(H+ DK) + (F+ GK )'S+S(F+ GK) = 0 (1) 

8J/3K = 2D 'Q(H+ DK)X+ 2G^X = 0 • K= -(D'Qpf 'QH+ G|S) (2) 

9J/8_S = X = (F+GK)X+X(F+ GK) ' + FX+XF 1 + CNC 1 = _0 (3) 

A 

Inserting K from (2) into (1) 

(F-G(D l QD)~ 1 D'QH) , S+S(F- qD'QD^D'QH) - ■ SGfD' QD ^tfS 
+ H^Q-QP^'QD^p'QJH = 0 

A 

With y; given, this equation may be solved for ^S. Then K is formed 

A 

using (2) and X is found from (3). 
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APPENDIX C 


VARIATIONS OF THE CRITERION FUNCTIONAL 


Consider a small change in X““ call it Av. Demand that 
Eqs . 1-3 of Appendix B be satisfied at all _y, i.e. , at y; an d 
This means that the first and second order variations of the above 
equations must be zero for all Ay, yielding variational equations in 
S, K, and X. 

First variations of Eqs. 1 and 2, Appendix B: 


(F + GK) 'S.. + S^.(F+ GK) + FjdS + SF^ + (Hj. + D^K) '0(H + DK) 

(1) 

+ {H+DKI'Q^.+p^K) = 0 

D l QD K li + D|.Q(H + DK) + D'pff^ + Dj^K) + G'S^. - 0 (2) 


Equation 2, Appendix C has been used to eliminate _X^. and 
. terms which are not shown above. 

l 

The second variation of Eq. 1, Appendix B is: 


(F+S® 'S 2 . .+S 2 . .(F+GK) + (Fjj+GKjj)' Sj. tSytFj.+ GK^ 

+ (F lj+ G5,j)'S u+ Sj.tFj.+GKj.) + (H li+ D 1 .K + DK li ).Q(H 1 . + D 1 .K + pK 1 .) 

+ (H lj +Sj j S+2K lj ) 'S (Bji+Dj.K+DKj.) + lF 21j+ GK 2 . .) '£ 

+(H 2ij +D 2i .K + DK 2 .. + D 1 ^ 1 .+D ]j K u ) '8(H+DK) + S(F 2iJ+ GK 2i .) 

HH+DK) 'Q(H 2ij +fi 2ij K+DK 2 . . + Dj.K,. + Djjy =0 (3 ) 

Using the preceding equations , this can be manipulated into the 
following form: 
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2ij + — 2ij^— + — — ^ + — 2ij-*— ^2ij + — li%j + — lj— li + — lj— id— li— lj 


+ (H 2 . .+D 2 ..K) 'Q (H+ DK> + (H+DK) ' Q (H 2i . + D 2i .K) 

+ (H li +P 1 .K)'Q(H lj +P lj K)+ (HJ.+ D lj K)'Q(H u + P li K) 

- Kj.p'QD^ - Kj j P ! QDK li = 0 (4) 

In the general case, one takes two variations of Eqs. 1, 2, and 

/S 

3 from Appendix B, yielding six equations in S, K, X, S^, K^,, X^> 
a 2 

S,.., Kt. and X,... This is a total of 3(1 + n + n ) equations which 
-2ij’ — 2ij -2ij y y 

may be solved for all of the variables above. Closer investigation re- 
veals that the solution procedure for the variational sets of equations 

A 

is the same as for the S, K, and X equations, except that the 
variational equations are of state covariance rather than Riccati form. 
The n^. + 2n^ + 3 equations given above (Appendix B and (1), (2) and 
(4)), however, are all that are needed in the remaining derivations. 



APPENDIX D 


CONDENSATION OF THE SECOND VARIATION EQUATIONS 


In order to carry out the Newton -Raphs on search it is necessary 
to compute the first and second variations of the criterion; this ap- 
pendix contains a series of manipulations which reduces the number 
of equations needed to compute these variations from a potential 
3(1 + n^ + n^.) to 3 + 2n^, minimizing computational effort. 

Using Eq . I of Appendix B, rewrite Eq. 1 5 of Section II as 
follows: 

J(i) = TR[S(F X + XF + CNC l ) + H|QHX ] (1) 

Because of the form of X (Section II, Eq. 13), _S is the only term of 
the criterion having a y dependence. Hence the first variations of 
J{y) are: 

J u = TR[S^(FX + |f' + CNC’)] = -2TR[S 1 .(F+GK)X] (2) 

byEq. 1, Appendix B. Similarly, the second variations maybe 
written 

J zi . = -2TR[S 2 ..(F+GK)X] = TR[{-S 2 . (F+GK)-(F+GK)'S 2ij )X] (3) 

But by Eq, 4 of Appendix C, these are 


J 2y = 2TR{[ (SF 2 .. + S iJ F li+ S 1 .F. j+ (H+DK) D 2ij K) 

+ (Hy + DyK) ’OfHj.+Dj.K) -KjjD'QDKj.lJX} (4) 

A 

Thus computation of the variations requires solution of the S_, K, X, 

and Kj. equations, or one Riccati equation, n^+1 state covariance 
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equations and Hy+1 direct evaluations . (Note: K^. and K. . may be 

eliminated from the J^ij ec luations, but no essential simplification 


results . ) 



APPENDIX E 


EQUATIONS OF RIGID BODY MOTION AS 
LINEARIZED EULER’S EQUATIONS 

The following derivation is a synthesis of Landau and Lifshitz' 
treatment of "Motion of a Rigid Body," and Etkin's "General Equations 
of Unsteady Motion, "(Refs . 1 and 2). Mechanics provides the basis 
for the notational conventions (with 1,2,3 corresponding to Cartesian 
coordinates x, y, z and angular momenta being referred to the axes 
about which rotation occurs), while Dynamics of Flight is used to 
phrase the problem in aerodynamicists 1 terminology. The derivation 
attempts to provide motivation for the linearized equations of motion 
on the basis of first principles well-known in classical physics and 
to put the flexure problem in proper perspective to the over-all control 
problem. Since this appendix is intended only as a link between physics 
and aerodynamics for the control engineer, details in the derivation of 
aerodynamic and propulsion force terms are not given. 

It is well known that the motion of a rigid body can be described 
by (1) the motion of its center of gravity with respect to a fixed co- 
ordinate system (e.g. , fixed to the earth or fixed with respect to a 
mean wind), and (2) the change in orientation of the body with respect 
to its center of mass. Newton's second law governs the three trans- 
lational and three rotational degrees of freedom: 

Ref 7 1 . Etkin, Bernard, Dynamics of Flight , John Wiley and Sons . , 
Inc., New York, N. Y. , 1959, Ch. 4. 

Ref. 2. Landau, L. D., and Lifshitz, E. M. , Mechanics , Addison- 
W esley Publishing Co. , Inc., Reading, Mass., I960, 
Pergamon Press, Ch. VI. 
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(1) 

(2) 


d'P/at' = F 
d 'M/ dt 1 = K 

P = linear momentum of the aircraft (time derivative 
with respect to fixed axis system) 

- mV, where m is the mass of the body and V 
is the velocity of the center of mass 
F = sum of forces on the body 

M = angular momentum of the body (time derivative 
with respect to fixed axis system) 

= I • £2 , where I is the inertia tensor and £2 
is the angular velocity of the body about its 
center of mass. (Note: I becomes diagonal 

when the principal axes of the body are used) . 

K = sum of torques on the body 

Because M bears a simple relation to the angular momentum with re- 
spect to the center of mass of the body, it is desirable to express the 
time derivatives above with respect to the body co-ordinates rather 
than the fixed co-ordinates; letting (d/dt) denote differentiation with 
respect to the body co-ordinates, one writes for any vector A, 
d’A/dt' - dA/dt + £2x A (Note: this convention does not follow Mechanics) . 
Thus (1) and (2) become: 

dP/dt + flxP = F (3) 

dM/dt tQxM = K (4) 

If the body axes are chosen as principal axes, these six equations may 
be written out in component form as follows: 

dVj/dt + J2 ? V 3 - S2 3 V 2 = Fj/m 


(5) 
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dV 2 /dt + 

= F 2/ m 

(6) 

dV 3 /dt+S2 L Y 2 -J2 2 V 1 

= F 3 / m 

(7) 

^d^/dth (I 3 -I 2 )f2 2 f2 3 

11 

J* 

(8) 

I 2 dS2 2 /dt + (I 1 -I 3 )fi 3 ^ 1 

= K 2 

(9) 

I 3 dft 3 /dt + 

= k 3 

(10) 


These are Euler's equations. (Note I. = (I - ) jj) 

The next problem is to choose a convention for expressing the 
rotation of the body co-ordinate system (x,y, z) with respect to the fixed 
co-ordinate system (X, Y, Z) and to express the components of the 


angular velocity of the body (£ 2 ^, £ 2 2 , ^ 3 ) terms of the change in 
rotational position. Euler's angles define the convention to be used. 



Fig. 1 Euler Angles 
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To rotate the xyz axes from a position of alignment with the XYZ sys- 
tem to a position of alignment with the principal vehicle axes, one 
carries out three rotations, in the following order: 

(1) A rotation tj) about the Z axis 

(2) A rotation 6 about the position of the y axis re- 
sulting from rotation (1), and then 

(3) A rotation <j> about the position of the x axis re- 
sulting from rotation (2). 

For clarity of presentation, the final position of the aircraft above was 
chosen in such a way that 9 and cj> as shown happen to be negative. 
Aerodynami cists usually define the z axis "down"; note that the 

signs of the angles follow the right-hand rule. 

• . * 

By considering the directions of ty, 9 , and <{> relative to the final 
position of the xyz system along the principal axes of the vehicle, one 
obtains the following expression for the components of Q : 

S2 ^ = (j) - ^ sin 6 

= 0 cos <{> + ip cos 8 sin <j> 

* • 

^2 = $ cos 8 cos 4> - 0 sin <j> 

If the principal axes were to be used in a linearized analysis , one 
would simply proceed to linearize Euler's equations as given above; 
this is sometimes done. More often, however, one uses a set of 
"stability axes." First, a condition of steady forward flight is speci- 
fied about which to linearize; the flight condition cannot be chosen 
arbitrarily, but corresponds to a detailed b alance of the aerodynamic 
and thrusting forces on the vehicle for a given "mean wind." In order 
to achieve a constant altitude at a specified speed, the nose of the 
aircraft must be tilted slightly upward (or downward) by an angle, say 
a Q (Etkin uses e) . Hence, in a given flight condition the velocity 
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vector of the center of mass and the principal axes don't coincide; 
one is thus motivated to use a set of "stability axes" aligned with 
the forward velocity vector of the aircraft, rather than the principal 
axes . 


In the following derivations the most commonly encountered 

flight condition is used as an example, viz., steady forward speed, 

Y^, constant altitude and small a Q , with no mean wind. With stability 

axes, the inertia tensor has nonzero components: 

2 2 

I, , = I. cos a + I, sin a — L 

11 1 o 3 o 1 

*22 = *2 

.2 2 

1, 0 = L sin a +1, cos a — I_ 

33 1 o 3 o 3 

and I., = I-, =--y{L-IJ sm2a - a (L.-IJ 

13 31 Z 1 y o o 3 1' 

where single subscripts refer to the principal moments of inertia. The 

reference values of V^, , ij), d 9 and <j> are zero. 

Carrying out the linearization of Eqs . 5-10 (as modified by the 

presence of L^), using small letters for the perturbed values, one 

obtains : 


dv^/dt = f-^/m 

(ID 

dv^/dt + ip = f^/m 

(12) 

dv,/dt - V, 0 = £,/ 
y 1 y m 

(13) 

% 4> + a 0 ( I 2 _I i)^ = k x 

(14) 

l J = k 2 

(15) 

a o^3 _I 1^ + h $ = k 3 

(16) 


The basic equations of motion for a general linearized rigid body 
problem are thus six in number: horizontal, lateral, and vertical 
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acceleration equations; and roll (4>) , pitch(0), and yaw (ip) acceleration 
equations. The role of the aerodynamicist in such linearized control 
problems is to specify the linearized forces, £., and torques, 1c., 
i = l,2,3 which excite the vehicle. The forces fall into several cate- 
gories : 

(1) aerodynamic forces 

(2) gravitational forces 

(3) thrust and propulsion forces 

(4) gust forces 

(5) control forces 

Gravitational forces (usually negligable in the Incremental model) 
and aerodynamic forces are usually specified in terms of forward 
velocity, v^, angle of attack (a=(v^ and sideslip angle 
( (3— (v^hVj $/V^) , and the roll, pitch and yaw angles and their derivatives. 
Thrust and propulsion forces are usually related to throttle position. 
Gust forces depend on side, head-on and vertical gust strengths, which 
are usually stochastic inputs. Control forces are generated as func- 
tions of the aerodynamic surface deflections. 

Generally, the control engineer is given the coefficients of the 
linearized equations of motion and hence need not be concerned with the 
detailed physics of their derivation; the control surface location prob- 
lem, however, poses an exception to this rule. It is necessary here to 
estimate all types of forces which depend on control surface location. 
For the rigid body equations, the only effect of surface location is to 
change the torque arms of the control force contributions to kj, k^ 
and k^ . The other effects of surface placement arise from the flexible 


nature of the aircraft. 
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The relation o£ the rigid body equations of motion to the flexure 
equations is not as clear-cut in practice as it is in theory. Certainly 
the deflection of a mode gives rise to aerodynamic force and torque 
terms in Eqs . 11-16. The issue is whether, when and how to include 
inertial coupling between rigid body motions and motions of the flexi- 
ble vehicle. Most theoretical treatments claim there is no inertial 
coupling of rigid body and flexure motions; this is because the normal 
modes are defined as the responses with distributed force and moment 
terms of the flexure equations set to zero. Such modes, however, 
are non-physical. Because mode excitation must give rise to aero- 
dynamic forces, which in turn feed back to excite the modes, it is 
unreasonable to set such terms to zero; rather, only applied forces 
(such as those arising from control surface deflections) must be set to 
zero. This procedure yields the true mode shapes for a given flight 
condition, but there is still the problem of normalization of the mode 
shapes. Fortunately, Stilley and Pollack (Ref. 17) have shown that 
(at least for vertical deflection modes) it is always possible to define 
mode shapes in such a way that inertial coupling to rigid body motion 
can be eliminated. Because normalization conventions vary widely, 
it is necessary to decide for each application whether or not the rigid 
body equations of motion do or do not contain acceleration due to 

flexure modes. The computer programs assume that flexure acceler- 

• . « « 

ations are included by definition in the model variables v^, <j> and 9 . 

Comment: Tbe equations of motion used in the C-5A model cor- 

respond to Eqs . 13 and 15. Since vertical gusts are the only inputs to 
the system and since the angle of attack is defined as above, these two 
equations decouple from the others, allowing one to treat the longi- 
tudinal axis separately. Decoupling of longitudinal and lateral equations 
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of motion is a general property of the equations of motion for the 

standard flight condition described above. For the C-5A equations, the 

z equation of the original data was purged of flexure contributions, 

* ♦ 

but 6 was not. 



APPENDIX F 


ORTHOGONALIZAT ION OF FLEXURE MODES 


Mode Acceleration Equations: 
oo 


oo 

Y, + X (Eie i" ),, Tli = F z {Y>t) 


m 


i=l 


oo 


i=l 

oo 


+ 2_^{GJ«j) i , ) , Ti i = t { y , t) 
i=l i=l 

The Normal Modes <j>. satisfy the equations: 


’or 


m(g i -S(j>.){or+ (Elg. ) = 0 

+ (GJ4>!)' = 0 


{ ’ =d/<l y > 


I . I 


= U?ms5. - (GJ<j>!) 


co?m£. = a/ms(|>. - (EI^!')" 


/end /end /end 

f / 1(4).^.)^ = JJ ms 5.<j).dy -J ( 


So wf* / I(cj) 
0 

y end 


(GJ <j>!)(|)!dy 


y end 


end 


(GJcjy^.dy 


w j/ I ( t l > i < l ) j) d y = w j/ msg^-dy-J* 

0 0 0 


■end "end 

**• (co. -an ) j !{<|>.<|>.)dy =J S-^msdy 


"end 


-/ [ (GJ$l )'<j>. -(GJcjy '<!,.] dy 


(1) 

(2) 

(3) 

(4) 

(5) 

( 6 ) 
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But as a boundary condition 


GJ *j (y end> = 0 ( freee ^ da t y end ) 
- = 0 (fixed end at y=0) 


On integration by parts, the last term is: 


[ (GJ«t»!)4>- - (GJ<j>'}^] 0 


end 


y end 


3 

^end 


for all j 


= 0 


•* {w 2 -w 2 ) J If^.Jdy = J (wj £.<|>. ~ w^j4>.)ms dy 


0 


(7) 


By precisely analogous manipulations on the Eq * 3, one establishes 


that: 



y end 


w j) f m ^i^j d y =/ (, -°i^i^j 


a). £ .<b.)ms dy 
J r 1 


(8) 


Adding, one obtains: 


y end 


M ij =f - ms(4>.g ; .+ cj>.^.)+ I<t>.4).] dy = 0 i/j 

0 

This is the orthogonality condition in the case of structurally coupled 
bending and torsional modes. 


When i=j 


y end 


M 


-jj = M j =/ - 2ms<j>.^. + I<{£] dy / 0 


3 J J 


Since I = Iqq + ms 


y end 


y end 


M. = 
J 


/ m[£j " s<|)j] 2 dy+J* I CG <f>j dy 
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Now, multiplying (1) by and (2) by <{k, integrating over the length 
and adding, one finds : 


f - ^(^^+^4.) 4- I<j>.({>pdy 


i 0 


'end f ena 

+T) i J [ (Elg.")"^.!- ( GJ 4> !) cjjj] dy = j" (F z (y,t)^.+ T(y,t)<t> j .) dy 


Referring to Eqs. 5 and 6, there results 


ena 

+ w? r|.) = f (F z (y,t)^ + T(y,t)(j)j)dy 


Or since M. . = M.6.. 

ij J iJ 


*end 

M M (t) + M coJr| (t) = J (F (y,t)£.(y) + T(y, t )<{>.(y)dy 
J J J J J * Z 3 J 

0 


as claimed in the text. 



APPENDIX G 


DESCRIPTIONS OF THE MAJOR SUBPROGRAMS 


G.l SUBROUTINE POTTER 

This subroutine solves the steady-state Riccati equation using an 
algebraic algorithm discovered by J. E. Potter. In order to solve the 
equation F^X + X F^ - XQjX + - _0 one first forms Potter's matrix 

ZP, where 

r* 

Q, 


zp = 




— i -?i 


If the equation is of n^ 1 order, ZP is a Zn x 2n matrix. Next the 
eigenvalues and the eigenvectors corresponding to eigenvalues with 
positive real parts are found; the roots of ZP can be shown to be the 
poles of the closed loop control system and their negatives, so there 
will be exactly n roots with positive real parts, assuming no zero 
roots. The n eigenvectors corresponding to these roots (real and 
imaginary parts of eigenvectors for complex roots) are then arranged 
column -wise and the resulting n x 2n matrix is partitioned, so that: 

A , where A is the Jordan real form for 

for the eigenvalues of ZP with positive 


X 

Qz“ 


>r 

II 

bf 

Pi 



2z m 




real parts . 


-1 


Then the solution may be written X = « Proof: 
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Right -multiply (1) by Y^ 1 ; left -multiply (2) by Y^ * and right- 
multi ply by Y^ , gives: 


- 1 ) + Q 2 

= ’pxf 1 

(3) 


= IiMl 1 

(4) 


Subtracting (4) from (3) and letting X = Y^Y^, 

FjX + XF^ - XQ]2£ + Q z - 0, as desired. 

Potter has shown that X will be the positive definite solution of the 
problem, providing is positive definite. 

The program forms _ZP and then uses an eigenvalue -eigenvector 
routine, ALLMAT,* which utilizes the double QR iteration of 
Francis 4 to find Y^ and Y^ . X is found directly by solving the 
system Y ' X = Y using the Gauss -Jordan method. 

The primary appeals of Potter's method as compared to iterative 
schemes are (1) speed and (2) accuracy. The savings in speed goes 
roughly as the system order {e.g., a 10 to 1 saving for a tenth order 
system), and the accuracy appears to be considerably better due to 
far more moderate error propagation. Note also that no free param- 
eters (such as step size) are involved. To the author's knowledge 
this method is the most efficient technique available for solution of the 
steady-state Riccati equation; the computer program is available on re - 
que st . 

G.2 SUBROUTINE STCOV 

This subroutine solves the state covariance equation AX+XA'+S: 

0. for X, where A, X, and _S are N x N and S and X are symmetric. 

jJU 

Courtesy of Union Carbide Corp. (Nuclear Division), Oak Ridge, 

Tenn. (Authors: John Rinzel and R. E. Funderlic). 
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The condition for convergence of the iterative algorithm (i.e., for 
existence of a steady- state solution to the associated differential 
equation) is that the eigenvalues of A have negative real parts; the 
program assumes this to be true. Certain convergence checks are 
provided, however. First, one may specify the maximum numbe r of 
iterations via the parameter ITMAX appearing in the subroutine CALL 
(if the first time increment is ~ l/3 0 of the system settling time, 
ITMAX - 15 has been found sufficient for good convergence). Secondly, 
one may demand that the fractional change of all diagonal elements of 
X at the final iteration be below the constant CTMAX appearing (cur- 
rently as . 001) in the program listing. The vector TRX is used to 
store the values of diagonal elements from the previous iteration and 
may be printed out as an exact check on convergence if desired. The 
remaining parameter, FR, is set to be the fraction of system settling 
time taken for the first iteration; settling time is approximated as 
TSET = l/MAX(A. ). Stein* has shown that the error in the initial 
iteration is always propagated as a constant in the answer and will 
never lead to a "fictitious" instability in the algorithm itself. The in- 
put matrices A and _S are not preserved. Computation time is well 
under a second for a tenth order example (using the IBM 360-40 com- 
puter) . 


* 

Honeywell Interoffice Correspondence (MR10033), G. Stein to 
V. Levapi, "A Note on High-Speed Computation of Steady State So- 
lutions to Linear Differential Equations, " Systems, and Research 
Center, July 14, 1967. 
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The method is an iterative algorithm for integrating the dif- 
ferential equation associated with the system above; the computational 
speed of the algorithm is due to the fact that step size is doubled with 

each iteration. The initialization and generation sequence are as follows: 

-A- A t 

Define _Z(1) = e— ~ [I + AAt), where At is chosen to be some 
fraction of the settling time of the system, approximated by 
(l/max (A..))). 

M 3 

At 

Define X{1) = j Z(t-T)S_Z'(t-r)dt 

0 

~ (At/2)((I + AAt) + i_I+ A'At)) 

Then iterate using the algorithm: 

Z(2i) = (Z(i)) 2 
X (2i) = X(i) + Z(i)X(i)Z'(i) 
i = 1, 2, 4, 8, 16, . .. 

The process ends when the desired convergence (via one of the cri- 
teria above) is attained. 

G.3 SUBROUTINES GD, GPDP AND GPPDPP 

Subroutine GD directly computes all elements of G o (_y) and D(y) 
given and the parameters of the equations in Sections III. A. 3 and 
in. A. 4. The basic pattern is to identify the equation, identify the matrix 
elements which it contributes, and then proceed to compute these ele- 
ments. The forcing terms of the mode acceleration equations are com- 
puted first and then added into the rigid body equations (those which are 
not ortho gonalized to the flexure modes) as the other ^-dependent terms 
are computed. As shown in Section III, the modal forcing functions and 
the rigid body forcing terms can then be used to compute the terms of the 
response matrix D(y). 
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Definition of variables: 

TSEX = conservative approximation of system settling time 

CTMAX = maximum allowable fractional error between estimates for satisfactory 
convergence 

DT = At for first iteration 

FR = fraction of settling time for first iteration 

TIME = time to which equation has been integrated 

NIT » current number of iterations 

Z - transition matrix corresponding to the system X = AX (Z is not a matrix in 
the program itself but is used as an identifier in the flow chart to show 
which storage locations contain che current transition matrix^ this is the 
same as the Z used in the thesis proposal write-up). 

X = current value of solution matrix 

TUX = trace of previous solution matrix 

CT = fractional change of diagonal element for current iteration 

NEZ - number of diagonal elements converged so that CT < CTMAX 

ITMAX = maximum number of Iterations (run time) allowed for convergence 

N = system order 


Fig- 1 Subroutine STCOV (Contd.) 
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NOTE: The left-hand side of an equation represents the storage locations in 

which the right-hand side is stored. For example, A = Z{l) means that the storage 
locations originally assigned to the A matrix contain the elements of the initial 
transition matrix, Z(l)- Considerable shuffling occurs in order to save storage 
space. 
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The variables used by the subroutine may be broken up as follows: 
(1) G q and D matrices, their dimensions and constant terms; (2) de- 
scriptors of state, response, and control vectors; (3) variables speci- 
fying force producer position; (4) aircraft geometry parameters; (5) 
mode shape variables; and (6) lift curve parameters. A list of vari- 
ables follows. 

(1) G o andD matrices, their dimensions, and constant terms; 


( 2 ) 


G(NX, NU) = 
DY(NR,NU) = 
NX = 
NR = 
NU = 

GSTOCK(NX, NU) = 


DSTOCK(NR, NU) = 


Descriptors of state, 


G (y) matrix 
_D(y ) matrix 
dimension of state vector 
dimension of response vector 
dimension of control vector 
constant terms for jy- independent 
part of G£y); positions of _y -dependent 
elements contain coefficients of these 
elements 

constant terms for y-independent part of 
D(y); positions of _y-dependent elements 
contain coefficients of these elements 
response and control vectors: 


IXPOS(NX) 


IXPOS(l) = 

position of z 

equation in state 


vector 


IXPOS(2) 

position of 6 

equation in 


state vector 


IXPOS(3) 

position of y 

equation in state 


vector 
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EX POS(4) = position of ^ equation' in 

state vector 

IXPOS(5) = position of <j> equation in state 

vector 

IXPOS(6) = position of x equation in state 

vector 

IXPOS(K+6), K= 1 , NI - position of in state 
vector; NI = number of flexure 
modes considered 

IUPOS(I, J) = position of force producer J of 

type I in the control vector; 1=1, 

NT; J=l, NMAX, where NT(=9) = 

number of force producer types and 

NMAX = maximum number of 

force producers of any type 
th 

IRPOS(I,J) = position of J response of type I 
in the response vector 
1=1- acceleration responses 
1 = 2- stress responses 
1 = 3- stress rate responses 
IRT = 3 = maximum value of I 
NRT = maximum number of responses of 
any type (maximum value of J} . 

YR(IRT=NRT) = response station location on 
airplane 

IRSEG(IRT=NRT) = axis upon which response 
station is located (i.e., y^ = 1 , 
y 2 = 2, y 3 = 3, y 4 - 4) 

NRI(I), 1=1, IKT = number of responses of type I 
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(3) Variables specifying force producer position 

LOC(I) ,1=1, NT - number of axis upon which I 1 * 1 type 
of force producer is located 
NFP(I) ,1 = 1, NT - number of surfaces of type I 
{see Section III for definition of types) 

YFP(I, J), 1=1, NT, J = 1, NMAX - location of J th 
force producer of type I, referred to axis system 
described in Section III 

AXANG(I), 1 = 1, NT - angle between center of pressure 
axis for control surface type I and fuselage 
AXX(I) ,1 = 1, NT - distance of root of center of 

pressure axis of 1^ type control surface fore of 
center of gravity 

(4) Aircraft geometry parameters 

YO(I), 1=1,4- start of I^ 1 aircraft (flexure model) axis 
YEND(I) , 1=1,4- end of 1^ aircraft (flexure model) axis 
YCG = position of center of gravity on y^ axis 
EAXX(I) , 1 = 1,4- AXX for I ^ segment of elastic 
axis system 

EAXANG(I) , 1=1,4- AXANG for I th segment of 
elastic axis system 

(5) Mode shape variables 

AK, AKK(I, J,K); 1=1, NI, J = 1,4, K = 1, NFIN - 
cosine coefficients in Fourier series for mode 
slope derivative of vertical (AK) and torsional 
(AKR) displacements for I** 1 mode, J** 1 axis, 
component in the series (NFIN = maximum number 
of terms in series) . 
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BK, BKK(I, J,K); 1=1, NI, J = 1,4, K = 1, NFIN - 
sine coefficients in Fourier series for vertical 
(BK) and torsional (BKK) mode slope derivatives 
for I*'* 1 mode. 

PERIOD(I), I = 1,4- period of fundamental in the 
the Fourier series for axis I. 

LAST (I) , 1 = 1,4- number of last term in Fourier 
series for axis I 

ZO(I, J-), ZPO(I,J), 1 = 1, NI, J = 1,4 - mode shape 

(ZO) and mode slope (ZPO) for vertical dis - 
th 

placement of I mode evaluated at YO of axis J. 
TO{I, J) , TPO(I, J) , 1 = 1, NI, J = 1,4 - mode shape 

(TO) and mode slope (TPO) for torsional displace- 
ment of mode evaluated at YO of axis J. 

til 

Z2M(I) , 1 = 1, NI - normalization constant for I cn mode 
generalized force (see Appendix I) 

(6) Lift Curve parameters 

CLAK(1, 4, NFIN) - coefficients of cosine terms in 

Fourier series for section lift curve (second index 
= axis number, third index = term number, first 
index always unity) 

CLBK( 1 , 4, NFIN) - coefficients of sine terms in 

Fourier series for section lift curve 

CLO, CLPO(l, J), J = 1,4 - value of section lift curve 

and lift curve slope at YO for axis J. 

CLC(I, J) I = 1, NT; J = 1, NMAX - coefficient by 

til 

which lift curve is multiplied for I type surface - 
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e*g. , the wing section lift curve will be multiplied 
by one constant for ailerons, another constant 
for spoilers, etc., the constants depending on 
size and position of the surface type. 

Working arrays in GD 

COS3, SXNS(NFIN) - used to store values of sine and 
cosine series terms 


RROW, RCOL, ZEV(NI) - used to store vector of NI 
modes evaluated at a given point 

CL(NT,NMAX) - used to store lift coefficients for each 
control surface 

Evaluation of the _y- dependent portion of a plant equation (say 
th 

the K equation) follows a pattern similar to the example 
below: 

NROW = IXPOS(K) = row of Cj£y) under consideration 

NCOL = IUPOS(I, J) = column of G (y) under con- 
sideration corresponds to position of control sur- 
face (I, J) in control vector, u. 

G(NROW, NCOLi) = F(YFP(I, J)), F being the ap- 
propriate function of 

NI 

G(NROW, NCOLi) = F(YFP) + C T (YFP) . S C.( YFP)£.{ YFP) 

L k=l J J 


where C.(YFP) is the appropriate constant (e.g. , 

J 

3 £ . 

dy~ ^Y cq ) f° r & > e t c >) f° r the modal acceleration 
terms . 


Generally there are three nested DO loops - one on I (type of con- 


trol surface), one on J (number of force producers of type I), and 
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one on K (for mode acceleration pickup in the rigid body equations). 

The DO loop on I usually skips terms to pick up either collective 

deflection contributions to pitch equations or differential deflection 

contributions to lateral equations. The subroutine EVAL is called 

several times to. evaluate 

C T , if these are functions of y, 

L. ■'k 


( 1 ) 

( 2 ) 

(3) 


§j(y k ) and J = 1 > NI 

9£- 

^j(yC G )» 8y ^CG^’ etc * usec * f° r tk e coefficients 


C.(y.) 

3 


Subroutine GPDP 

9G(i) SD(x) 

This subroutine evaluates the first partials — r and — - 

dy i 

where y.=YFP(ID, JD) is specified by ID and JD in the subroutine 

CALL. The derivatives are evaluated directly by differentiating the 

expressions given in Section III. A. 3: 

NI 


dz y /dy r (dC L /3y i )6 i + £ (^(y^HdX./dy^/ZZMfj) 
1 i=l 


90 /8y. = 
Y 1 


9 C 


i=l,3,5,7 

(1) 


L. 


9y 


~[ AXX(i)-y./TAN(AXANG(i))] -C T /TAN(AXANG(i) 

. 1 _L 1 . 


NI 

^-4 9g. 9X. 

+ / (y r.n> 


(2) 


L 9^ Z2M ^> 

j=l 2 


i = 1,3, 5,7 


9 Y /9y. = 
v' r 


0 


(3) 
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3 a 


3x /ay.^ -g~i- [ AXX 2 (i)+y 2 (l+COT 2 (AXANG(i)) 
y i y i 

-2y.AXX(i)COT(AXANG(i))] 


(4) 


+ C D C L. W- 


[ y.(l+ COT 2 { AX ANG{ i)) ) -AXX(i) COT( AXANG(i) ) ] 6 . 

“i "i-1 1 [AXX 2 (i)+y 2 (l+COT 2 (AXANG(i)))-2y.AXX(i)COT(AXANG{i))] 1 ^ 2 

i = 2,4, 6 


9 4-/^9 = C D COT{AXANG(9)) 

9 

3<p /3y. = C p _(y.3C L> (y^/Sy. + C L _ (y.))6. 


(5) 


1 

NI 


i-1 


i-1 


9 |,- 


j=l 


i = 2 , 4 , 6, 9 


34 > y / 3 y 8 


= C, 


9C L ^8 )/8y 8 (y VT +y 8 )l/2 + C ^ (y 8 )y 8 /(y VT +y 8 ) 


1/2 


j=l 


3X /ay. = C D 3G L /8y. 


^a£. ax. 

+ Z ^ y CG ^^ Z2M ^ 


X X 


- C D. 8C I.. /<*! 
1 1-1 


i = 1,3, 5,7 


k = 2, 4, 6, 8 


( 6 ) 


= 0 


The term 

ax. 


i = 9 


aC L./ 0y i^j( y i) + d ( y i> ( y i>]Si 

+ C L.( y i)t 9 ^j/ 8y i + d (y i )34> j /3y i + ad/ay.^^yp] 6. 


(7) 



-130- 


In the 


computer program d(y\) = DIS and 9d/8y. = DISP. Here g 


.th 


is the vertical deflection of the j mode (see above) and (j>. is the 

th 

torsional deflection of the j mode. The subroutine EVAL evaluates 

mode shapes mode slopes {9^j/9y., 8^/3 y.), and lift curve 

shape and slope (C^ , 3C^ /9y.) . 

i i 

The structure of this subroutine is basically the same as in GD 
and the variable names are unchanged. The major differences are 
that the first variation matrices GP = 9G^/9y^ and DP = 9D/9y^ have 
only one nonzero column corresponding to control surface i {note that 
a single index i has been used on y throughout this description, 
whereas in the program, y. = YFP(ID, JD) , is identified by two indices -- 
one for the type of surface and one for the number of surfaces of that 
type). In GPDP, the constants CLIFT = and CLIFTP = 9C r /9y. are 


L/ 


used in place of the matrix CL = which stored lift coefficients in 

i 

subroutine GD . As in GD, the constants of y- dependent terms 

(e.g., C D , Cp , etc.) are stored in GSTOCK and DSTOCK. A flow 
i i 

chart for the program is shown in Fig. 3 . 

The acceleration terms of 9D/3y. are proportional to: 

1 NI 

8a z / 8y i = 8 V 9y i " x r ' 8 ^ y / 8 > r i + y r 9 <f> y / a yi + X <0X j/ 8y ^j{ y r ) 

j=l 

( 8 ) 


The proportionality constant, again, is stored in GSTOCK. Note that the 

d6 /dy. terms apply to the pitch axis control surfaces (i=l,3,5,7) and 
y A 

the 9<j> /9y terms apply to lateral axis control surfaces (i=2 , 4, 6, 8, 9) . 
The moment arms x^ = XR and = YR(1, K) , and the £.{y^) are 
computed first; remaining terms have been computed for GP- 
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Fip. 3 Subroutine GPDP 




























1 oo 
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Subroutine GPPDPP 

This subroutine computes the second variations GPP - 8 G o /9^8y. 

j 

z 

and DPP - 9 D/dy.ay^ of C Q (y) and D(y) . Note that the equations 

given in the GD description contain no cross-product terms in y. 

and y., and that force producer locations are varied independently so 

j 

that y. and y. are unrelated unless i=j. Hence GPP = _0 and DPP ' =_0 
for i / j. For i = j, GPP and DPP each have one nonzero column (as 
in GPDP) which contains second partial derivatives of the rigid body 
acceleration equations, mode acele ration equations and response ac- 
celeration equations- These are written out as: 

NI 

hy/dyh a 2 c L /ayf 6.+ £ e k (y CG )(a 2 x k /9 yi 2 )/z2M(j) (D 

1 k=l 

i = l,3,5,7 

d Z 6 Jdy z =[a 2 C L /8y 2 [ AXX(i)-y.COT(AXANG(i))] 

- 2<9C T /9y.)COT(AXANG(i) )] 6. (2) 

J-i. 1 i- 

1 

ni di 

+ Z -8^ ( yCG )(s2 V 8 n>/ Z2M( j» 

k=l 


i = 1,3, 5,7 


a 2 y y /a 2 y . = o 

r: * 2 


(3) 


..,2, 2- 


a ^iyf = C D 9 /9yf[ AXX(i) +y7‘ ( 1+COT (AXANG(i))) 

y 1 i 1 


-2 y . AXX ( i) C OT (AX ANG( i) ) ] 1//2 

(y.(l+COT 2 (AXANG(i))-AXX(i)COT(AXANG(i))) 
+29C t/ %. 1 


(AXX(i) 2 +y?(l+COT 2 (AXANG(i)))-2y.AXX(i)COT(AXANG(i))^ 72 


( 4 ) 



-133- 


+C 


(1-fCOT (AXANGfi))) 

L i -1 ( AXX( i) 2 +y 2 (l+C OT 2 (AXANG(i) ) ) -2 y. AXX ( i) C OT (AXANG(i)) 1 / 2 

y 2 ( 1 - AXX ( i) C OT ( AX ANG ( i) ) + CO T 2 ( AX ANG ( i) ) 2 
(AXX( i) 2 +y 2 ( 1+COT 2 ( AX ANG(i))) -2y. AXX(i) COT(AXANG{i) )) 


^ 2 


i = 2, 4,6 

9 2 W3 2 y. = C p {2 • 8 C l (y.)/9y. + y^ 2 ^ (y^y 2 ) 6 i 

y i i-1 i-1 


NI 3? 

V 9g k 

L 9 yj 

*1 1 1 


(y CG )(9 2 X k /3y 2 )/Z2M(j) 


k=l 
i = 2,4,6, 9 
27 /„ 2 „ 2 


2 ., 2 


V 8y S = C P 0 (8 \ ly 8 l '' 8y 8 ),y VT +y 8 )l/2 

o 7 


+2{8C L 7 (y 8 ) / 8 y8 ) / (y VT +y 8 )1//2 


+ C L 7 {y 8 )[l / (y VT +y 8 )1/ ' 2 " y 8/ (y VT +y 8 )3//2] 


NI A£ 

V J!j 

Z , ay -, 

k=l 1 


(y CG ) 3 Z X.(y 8 )/9y 8 )/Z2M(j) 


9x /9y. = C~ 3 C, /3y. 
Y i D. L,.' } l 

li 

= c d. 92< ?l. 

1 1-1 


= o 

The mode forcing terms are 

> 2 „ /„ 2 r ^.2 ^ 2 


i = 1,3, 5, 7 
i = 2,4, 6, 8 
i = 9 


9 x j/ 9 yi = {9 C L /3y. [^(y.) + dfy.^fy.)] 

+ 29C L /9y.[9£ j /9y i + d(y.)9<j>./9 yi + Bd/By.fyy.)] 

+ C L (Y^t./dyl + 23d/9y.8<j>./ay i ]}5 i 

(NOTE: 9 2 d/9y 2 = 0) 


(5) 


( 6 ) 


( 7 ) 
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The response accelerations are 


3 2 a z /ay? = d 2 Zy/dy* - x r 8 2 0 y /9y^ + y r 8 0 y /9y. (8) 


NI 

+ ^a^/ay^y, 

j=i 


.) 


The organization of GPPDPP is precisely that of GPDP, and variables 

retain their definitions. The subroutine EVAL is also used to compute 

second derivatives of mode shapes and section lift curves . 

1 

G. 4 EVAL 
Equations Solved: 


N 


z (y) = Z Q + %y+ (a o /4)y 2 - a k (L/2Trk) 2 cos(2irky/L) 


( 1 ) 


k=l 

N 


L/2 irk) 2 s in ( 2 irky /L) 


k=l 


N 


My) = Z Q+{ a 0 / 2 )y + ^ a j c (^/ 2 'n'k)sin( 27 rky/L)-^b k {L/ 2 irk)cos( 2 irky/L) 


k=l 
N 


N 
> 1 
k=l 


( 2 ) 


N 


z M (y) = a- 0 /2 + a^cos(2irky/L) + b^sin(2irky/L) (3) 


k=l k=l 

Here a^_ and b^ are the sequences of Fourier coefficients in the 
Fourier series for z"(y), N is the length of the series, L is the 
period of the fundamental, and y is the point at which z is to be 
evaluated. The program simultaneously evaluates any number of 
modes (z), mode slopes, or slope derivatives at the desired evaluation 
point y. The point y maybe anywhere on the aircraft. 
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Note on computation method: 

The sine and cosine terms in the series are rapidly evaluated via 
the trigonometric identities 

sin(n+l)x = cos(x) cos(nx) - sin(x) sin(nx) (4) 

cos(n+l)x = cos (x)sin(nx) + sin(x) cos(nx) (5) 

The program may also be used for evaluation of other y-dependent param- 
eters such as section lift coefficient, bending stiffness, and torsional 
rigidity, 

G, 4 SUBROUTINE EVAL 

This compact, fast little program is the reason that mode shapes, 
slopes and slope derivatives are conveniently stored as coefficients of 
Fourier series. Given any point on the aircraft, EVAL. will return a 
vector of mode shapes, slopes, or slope derivatives for all flexure 
modes evaluated at that point. Hundreds of mode shapes maybe evalu- 
ated in a single second. The subroutine has a counterpart, SEFOUR 
which has been used in the Fourier de composition of mode shape data. 

The subroutine itself is amply annotated so as to be easily ‘’lifted 11 for 
other applications. A description of variables follows. 

AK = ’'tensor 11 of Fourier cosine coefficients for mode slope 

derivatives; first index is number of the mode, second 
index is segment of aircraft (1-wing, 2 -fuselage, 3 -vertical 
tail, 4-horizontal tail) and the third index is the term in the 
series (with as the first element). 

BK = ’’tensor 11 of Fourier sine coefficients. Series generally 

tend to have only the first cosine term and odd sine terms 
nonzero, though the program doesn’t require this. 

NI = total number of modes 

NFIN = maximum length of series for any segment of the aircraft 



136- 


PERIOD = 


LAST 


ZO 

ZPO 

NSEG 

YON 

YEY 

NFN 

ZEV 

COSS 


vector containing the period of the series for each aircraft 
segment (4 elements). The fundamental period should be 
about four times the length of the segment, 
vector containing number of last term of series for each 
aircraft segment. The number of terms for a given seg- 
ment should be such that the shortest period in the series 
is approximately the length of a panel used to compute the 
mode shapes. Inclusion of higher order terms introduces 
extraneous "noise" due to the lumped modeling of the air- 
craft. 

matrix of initial values of mode shapes for each segment. 
First index is mode number; second index is segment 
number (1-4). 

matrix of initial values of mode slopes for each segment. 

ZO and ZPO correspond to z(0) and z'(0) ofEqs. 1 and 2 
and are used to obtain constants of integration for the series, 
number of aircraft segment where modes are to be evalu- 
ated (1-4) 

value of coordinate variable where axis (corresponding to 
NSEG) begins 

value of coordinate variable where modes are to be evaluated. 

1 for mode shapes 

2 for mode slopes 

3 for mode slope derivatives 

vector of NI mode shapes, slopes, or slope derivatives, 

depending on NFN (output) 

vector of evaluated cosine terms (storage) 
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SINS 

PI 

PIOL 

YE 


= vector of evaluated sine terms (storage) 

= ir 

= 2-rr/PERIOD = "frequency" of fundamental 
= equivalent point at which series is evaluated (since Fourier 
coefficients are evaluated for = 0, it is necessary to 
shift the axis by YON in evaluating the series). 
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Fig. 4 Subroutine EVAL 


















APPENDIX H 


TRIAL MODEL COEFFICIENT DATA 

The parameters and coefficient data for the three-mode trial 
model of the C-5A described in Section IV is presented on Tables 1-4. 
Charts 1-4 show the mode shapes of the first, third and sixth modes 
used in the trial model. The dimensions of the model variables are 
as follows: 


Variable 

Dimens ion 

Name 

a = w 


( inches / sec . ) 

angle of attack 



(inches/sec,) * 

pitch rate 

T., 1=1 , 3, 

6 

(inches /sec . ) 

normal mode velocities 

t., i=l, 3 

,6 

(inches ) 

normal mode deflections 

p., i=l, 2, 

3,4,5 

(dimens ionles s) 

Kussner gust delays 

a 

g 


(ft./sec.) 

gust angle of incidence 

5 

sp 


(degrees) 

spoiler deflection 

6a^, 5a^ 


(degrees) 

aileron deflections 

6 

e 


(degrees) 

elevator deflections 

m^, m^ 


(g f s) 

accelerometer signals 

m 3 


(inches /sec. ) 

rate gyro signal 

s., i=l , 2, 3,4,5 

((1/(12) (32. 2) )lb. /in? 

stress responses 

a^, i=l, 2, 

3,4 

(g's) 

acceleration responses 

Spi=l,3, 

6 

(inches) 

mode shapes (vertical) 

i-l> 2, 

3,4 

(inches) 

orthogonal body coordinates 

<{>., i=l, 3, 

6 

(dime ns ionles s) 

mode shapes (torsional) 
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Sign conventions for trial model are positive for; 

(a) Trailing edge down for all control surfaces (Note: the 
"zero" deflection position for spoiler flaps corresponds to 
the flaps deflected 30°, or half their full range, to permit a 
valid linearization) . 

(b) Wing tips down and aft for bending modes . 

(c) Trailing edge of wings down for torsion modes . 

(d) Pitch nose -up for 6 . 

(e) De creasing altitude for z. 




Tublc 1 
Plant Model 


Pi 

n 

p l 

p 4 

p s 

m « 

*$,794# 

-#.7948 

-l. 73# 1 

-1.7311 

0 

-.35055 

-.39799 

-.39799 

-5, 6842 

-5.4842 

0 

2.1038 

.20.401 - 

20.401 

4.7342 

4.7342 

0 

.70142 

0 

0 

0 

0 

0 

0 

14.599 

14.599 

-5.1442 

-5. 1442 

0 

-1.9377 

0 

0 

0 

0 

0 

0 

um 

3.2118 

3.42 8 4 

3.42*9 

0 

.31323 

0 

0 

0 

0 

0 

0 

-2.2 

Q 

0 

0 

0 

.§578 

0 

24.34 

0 

0 

0 

14.042 

0 

0 

-4.44 

0 

0 

1.443 

0 

0 

0 

-53.24 

0 

33.544 

0 

0 

0 

0 

0 


0 

0 

0 

0 

i.o 

-1.154 



I 
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+ J Tj 

-.0014509 -.0036953 -. 0641 ? 

-.00040394 -.00002406 .0010934 

1,0 0 0 


-.0001 5241 -.41393 

-.00095952 -.17437 

0 0 

Ax 


Table 2 
Sensor Signals 

+ 6 \ P, P 2 

-.0033636 -.59212 -.039328 -.039321 

,00062301 , 33870 -.012414 -.012414 

0 0 0 0 


Pi 

-.00057125 


-.0022638 

0 


-.00057125 

-.002 2 638 
0 
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w 

o/t z 

T 1 

T 1 

t 3 

-.002259 

016123 

004345 

.61545 

-.040193 

-.00494 

.010717 

.002208 

.90554 

-.056853 

-.009193 

-.011387 

.000546 

.020272 

-.007464 

.011624 

020647 

-.003803 

12533 

-.003 639 

-.032271 

- 056583 

.015456 

020311 

- 022388 

-.060422 

-.02 609 

. 000993 

048292 

-.033856 

-.05541 

-.03890 

.010905 

22427 

-.027222 

-.05369 

-.05469 

.01861 

34491 

-.03 8589 

- 063589 

- 09943 

03 5119 

55015 

- 072403 


Tabic 3 

Response Equations 


t 3 

t 6 

t 6 

p l 

P 2 

-1 9026 

-.018 744 

-2 6249 

- 09278 

-i 09278 

3.6342 

-.005142 

4 2301 

- 12055 

- 12055 

-1.2465 

005044 

1 7633 

-.88974 

-.88974 

-l 3265 

.052461 

4.2468 

-.02778 

-.02778 

22195 

.071748 

-1.1886 

.001542 

.001542 

-5 0565 

.028502 

8 8001 

-.41993 

-.41993 

- 65968 

014487 

-7.0216 

-.45400 

-.45400 

- 5 6 53 6 

028181 

-9 3949 

- 36562 

-.36562 

-2 6920 

.097 791 

. 1 7274 

-.19103 

-.19103 




Stresses m l/(12)(32 2) lb/m? 

Accelerations in g f s 

Control Displacements in Degrees 


P 3 

p 4 

P 5 

a 

8 


6 . 
0 

S 

S 

6 sp 

10836 

.10836 

0 

-.046017' 


71742 

. 15419 

15419 

-.025543' 

085468 

.085468 

0 

-.024255 


.57011 

.25621 

.25621 

14172 

-.044539 

-.04 4 53 9 

0 

.82391 


-3.5580 

-1.2710 

-1 2710 

-1 2953 

.14917 

. 14917 

0 

.0)7830 


.014234 

-.016017 

- 016017 

-.068960 

087109 

,087109 

0 

- 002485 


-1.1601 

-.000249 

(.092033) 

- 000249 
(.092033) 

.020248 

-.14794 

-.14794 

0 

-.157222 


-.96555 

-.63 825 
(.25722) 

-.63 825 
(25722) 

-.19986 

-.2053 

-.2053 

0 

-.005648 


-1.3576 

- 21854 
(.23 887) 

- 21854 
( 23887) 

-.37799 

-.31980 

-.31980 

0 

041112 


-2. U 92 

-.14305 
( 020241) 

- 14305 
( 020241) 

-.31141 

-.64786 

- 64786 

0 

.075801 


-4 2917 

-.072513 

- 072513 

- 21035 


+ Dc 
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Table 4 

Auxiliary Matrices 


ZO = 


.0973 

-.0223 

-.1663 

-.2865 

.0349 

.3051 

. 1687 

.3215 

.2012 

-.3578 

-.2575 

- . 6263 


ZPO = 


IRPOS 


NRI 


0, 

0, 

0 . 


.7615E-04 - . 1414E-03 -.6160E-04 

.3 132E-03 .4467E-03 .2696E-03 

. 6342E-03 - . 1612E-02 -.1097E-02 


IRSEG = 


1 

2 

1 

0 


7 

2 

0 

2 

1 


8 

3 

0 

2 

2 


9 

4 

0 

2 

2 


0 

5 

0_ 

0 

4 


0 0 


0 0 



"101“ 


’4959.2 " 

LAST = 

97 

PERIOD = 

9910.8 


41 


1497. 84 


3 7_ 


JL460 . 8 _ 


4 

5 
0 


YR = 


395. 

1369. 

1804. 

2406. 

120. 

746. 

1106. 

1804. 

0. 

0. 

0. 

0. 


0 

20 

0 
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Table 4 (Contd. ) 
Auxiliary Matrices 




-146- 


Table 4 (Contd.) 
Auxiliary Matrices 



60. 0 


”l299.8 "" 

YO = 

277. 5 

YEND = 

2755.2 


405. 57 


780. 03 


20.3 

_ _ 


385.5 

L, J 
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(IN.) 



Chart 4 


Vertical Deflection of Horizontal Tail 



APPENDIX I 


COMMENTS ON FORCE MATRIX PARAMETERS 

This appendix contains miscellaneous notes on some of the more 
esoteric details of the modelling process, 

1. Note on calculation of AXX, AXANG: The axis defined by 
AXX(I) and AXANG(I) corresponds neither to the hinge line nor the 
edge of control surface I. The axis should be calculated to correspond 
to the aerodynamic ’'center of pressure” due to the forces created by 
deflection of the surface. The book Foundations of Aerodynamics by 
Keuthe and Schetzer contains a derivation of center of pressure location 
for a flapped airfoil section (Sec . 5 . 5-5 . 8K If E is the ratio of flap 
chord to local chord (c(y)), then the center of pressure is a distance 
x - c(y) . K from the leading edge, where 

TT-S 

This has been derived from the expressions given in the above source. 
As E ranges from 0 to 1 the center of pressure starts at the 
quarter chord, moves back to about 35 percent of chord and returns 
to the quarter chord. 



K 


Fig, 1 Center of Pressure as Percent of Local Chord 
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The above source also gives the theoretical derivation of lift and moment 
coefficients which can be used in the linearized force equations if 
coefficients cannot be derived from the equations of motion for fixed 
locations. The magnitude of the vertical force exerted at the center of 
pressure due to a flap is 

L. = C j jqs(a-a LQ ), where 

= [ 2{TT-sin -1 (EN/i7E)) + 4 En/L-E] 6 f 

6^ - flap deflection (radians) 
g = y p - dynamic pressure 

s = c(y) = flap area 

J?£ = flap length (span-wise) 

c(y) = local chord 
a = local angle of attack 
a ~ = zero lift angle of attack 

IjU 

2 . A few comments on the simplifications in the mode shape are 
in order. First there is the distinction between free aircraft modes 
shapes (|^) and the manufacturer's mode shape data for a specified 
flight condition. For a given flight condition there are distributed, 
velocity-dependent aerodynamic lift and damping terms on the right of 
Eqs . 1 and 2, Section III. A. 2. These terms arise because of a change 
in local angle of attack and hence air flow (lift) pattern over the wing 
when a mode is excited. These are.. customarily brought to the left 
side (being constants of each flight condition) before mode shape data 
is generated. With aerodynamic effects included, the "modes" are 
coupled and the "natural" frequencies are changed slightly due to 
damping terms . 
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Aerodynamic terms are augmented by structural damping and 
coupling terms between torsion and bending modes (these cause further 
aerodynamic terms). Structural damping is usually approximated for 
lack of sufficient data by a constant on the order of . 01 to .03. The 
structural coupling between bending and torsion modes depends on the 
specifics of basic structure of swept -wing aircraft. These structural 
effects are included in the equations for gi^en flight conditions and for 
the free aircraft. 

The net effect of aerodynamic and structural terms is to change 
the raw mode shape data given to the engineer and to add damping 
terms to the normal co-ordinate equations for the T^(t). Since these 
effects enter as data inputs (1 . e . , entries of the F matrix) rather 
than as computational modifications, they need be of no further concern. 
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